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Abstract 

Large,  lightweight  space  structures  of  the  future  will  require  state-of-the-art  vibration 
suppression  systems.  To  design  such  a  system,  it  is  necessary  to  have  a  mathematical  model 
that  adequately  describes  the  motion  of  the  system.  Coulomb  damping  compensation  was 
utilized  to  remove  a  known  nonlinearity  from  the  frequency  response  functions  obtained  from 
the  Passive  and  Active  Control  of  Space  Structures  (PACOSS)  Dynamic  Test  Article  (DTA). 
Via  the  Eigensystem  Realization  Algorithm  (ERA),  the  frequency  response  functions  were 
then  used  to  generate  a  mathematical  representation  of  the  structure’s  dynamics.  Two  models 
were  obtained.  The  first  model,  without  friction  compensation,  acted  as  a  baseline.The  second 
was  obtained  using  friction  compensation.  Comparing  the  models,  it  was  determined  that 
active  friction  compensation  was  worthwhile  and  resulted  in  a  more  accurate  mathematical 
description  of  the  system. 
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A  STATE-SPACE  MODEL  OF  A  LARGE,  LIGHTLY  DAMPED 
SPACE  STRUCTURE  USING  THE  EIGENSYSTEM 
REALIZATION  ALGORITHM 


1.  Introduction 

Future  military  and  civilian  space  structures  will  be  large  in  size  due  to  design  require¬ 
ments;  however,  in  order  to  reduce  launch  costs,  these  structures  must  also  be  lightweight. 
These  characteristics,  when  combined,  lead  to  a  structure  likely  to  experience  vibrations. 
These  vibrations  could  hinder  the  system’s  performance  and,  if  left  uncontrolled,  may  lead  to 
failure  of  the  mission.  This  is  especially  true  of  those  systems  requiring  precision  pointing, 
precise  shape  control,  and  rapid  targeting  maneuvers  (3:11).  These  systems  will  require  state- 
of-the-art,  optimal  control  systems  in  order  to  remain  useful  over  their  entire  life  span.  To 
perform  an  optimal  control  design,  an  accurate  state  space  model  of  the  system  is  necessary. 

1 .1  Problem  and  General  Approach 

The  goal  of  this  research  is  to  develop  an  accurate  state  space  model  of  the  Passive  and 
Active  Control  of  Space  Structures  (PACOSS)  Dynamic  Test  Article  (DTA),  a  testbed  for  large 
flexible  space  structures.  The  DTA,  like  actual  large  space  structures,  possesses  high  modal 
density  at  low  frequencies  (3:15,17).  This  leads  to  a  need  for  a  system  identification  technique 
capable  of  identifying  multiple,  closely  spaced,  low  frequency  modes.  In  this  research  effort, 
the  Eigensystem  Realization  Algorithm  (ERA)  was  used  to  develop  the  state-space  model 
(2:1564).  This  identification  algorithm  works  best  with  low-noise  or  perfect  data. 

In  order  to  obtain  the  most  accurate  state-space  (linear)  model,  known  nonlinearities 
should  be  eliminated  from  the  system  to  the  maximum  extent  possible.  Coulomb  damping,  a 
nonlinear  phenomenon,  is  known  to  exist  in  the  reaction  mass  actuators.  A  secondary  goal  of 
this  research  is  to  calculate  the  Coulomb  damping  force  in  each  reaction  mass  actuator  and  to 
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remove  the  effects  of  Coulomb  damping  from  the  resulting  frequency  response  functions.  The 
Coulomb  damping  coefficients  will  be  calculated  using  free  decay  (14:125).  After  determining 
the  damping  coefficients,  a  nonlinear  compensator  will  be  implemented  to  remove  the  effects 
of  Coulomb  damping  from  the  frequency  response  functions,  which  will  be  used  to  generate 
the  state-space  representation  of  the  PACOSS  DTA. 

1 .2  The  PACOSS  Program 

Martin  Marietta  Astronautics  Group,  during  the  Passive  and  Active  Control  of  Space 
Structures  program,  studied  active  control  and  passive  damping  techniques  for  use  on  large 
space  structures.  The  ultimate  goal  of  the  PACOSS  program  was  to  demonstrate  the  benefits 
of  passive  damping  working  in  concert  with  active  vibration  control  on  large  space  systems. 

Based  upon  planned  and  conceptual  civilian  and  military  space  systems,  a  Represen¬ 
tative  System  Article  (RSA)  was  designed  and  configured  such  that  vibration  control  design 
methodology  could  be  better  understood.  The  RSA  was  not  designed  to  be  mission  specific; 
instead,  it  was  a  representation  of  several  missions  and  requirements  in  a  single  system.  The 
RSA  was  composed  of  seven  substructures:  a  ring  truss;  a  box  truss;  two  solar  panels;  an 
equipment  tripod;  and  antenna  dish;  and  an  equipment  platform.  To  help  ensure  dynamic  sim¬ 
ilarity  to  future  large  space  systems,  the  size  of  each  substructure  was  modeled  after  planned 
or  actual  hardware  for  future  spacecraft.  This  guaranteed  that  the  damping  mechanisms  devel¬ 
oped  under  the  PACOSS  program  were  applicable  to  many  of  the  future  large  space  structures 
(7:33,38). 

A  scaled  version  of  the  Representative  System  Article  was  constructed.  The  Dynamic 
Test  Article  was  designed  to  have  characteristics  similar  to  the  RSA,  making  it  representative 
of  future  space  systems.  Thus,  damping  techniques  developed  and  tested  on  the  DTA  could  be 
used  in  space  systems  with  confidence  (7:1).  The  test  article  included  the  same  substructures 
found  in  the  RSA  as  well  as  large  amounts  of  passive  damping. 

After  extensive  study  of  the  structure,  Martin  Marietta  concluded  that  a  combination 
of  passive  damping  and  active  control  mechanisms  would  benefit  future  military  and  civilian 
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space  systems.  The  Astronautics  Group  did  report,  however,  that  improvements  were  needed 
in  the  area  of  system  identification  of  damped  stmctures  with  dense  modal  spectra  (7:138,141). 

Upon  completion  of  the  PACOSS  contract,  the  DTA  and  its  supporting  hardware  were 
moved  to  the  Air  Force  Institute  of  Technology  (AFIT)  for  further  study.  After  delivery,  the 
test  article  was  modified  by  removing  the  box  tmss,  the  antenna,  the  equipment  platform, 
and  all  of  the  passive  damping  mechanisms.  Due  to  logistical  constraints,  the  box  truss  was 
replaced  by  three  mass  simulators.  Having  the  same  dimensions  as  the  originals,  two  new 
solar  arrays  were  constructed  without  any  passive  damping  treatments.  These  modifications 
were  made  to  the  testbed  to  approximate  a  lightly  damped,  large  space  structure  where  purely 
active  control  techniques  could  be  developed  and  investigated  for  vibration  suppression  (4:2). 
Active  control  techniques  are  being  investigated  because  of  the  added  launch  costs  that  are 
incurred  when  passive  damping  mechanisms  are  utilized.  This,  of  course,  is  due  to  the  added 
weight  of  the  passive  mechanisms. 

Since  its  arrival  at  the  Air  Force  Institute  of  Technology,  the  PACOSS  DTA  has  been  the 
object  of  three  research  efforts.  The  first  effort  at  AFIT  was  conducted  by  Capt.  Scott  George 
(4).  Capt.  George  conducted  an  extensive  modal  survey  of  the  test  article  using  finite  element 
analysis.  Lt.  Chad  Matheson’s  research  involved  applying  purely  active  feedback  control  to 
attenuate  vibrations  in  the  test  article  (6).  Lt.  Matheson’s  research  also  involved  the  creation 
of  the  first  mathematical  model  of  the  flexible  body’s  dynamics  since  its  modification.  Capt. 
Tony  Nash’s  research  effort  was  to  develop  a  state  space  model  for  the  test  article  using  time 
domain  techniques  (8). 
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//.  Experimental  Equipment 


The  PACOSS  Dynamic  Test  Article  and  its  supporting  hardware  now  reside  at  the  Air 
Force  Institute  of  Technology.  The  support  equipment  delivered  by  Martin  Marietta  includes 
the  following:  the  air-bearing  suspension  system,  the  motor  control  units,  the  Optima/3  array 
processor  and  a  Sun3/50  computer.  Additional  equipment  needed  to  continue  research  on  the 
lightly  damped  structure  was  obtained  by  the  Air  Force.  This  chapter  describes  the  DTA  and 
its  equipment  as  well  as  their  use  in  this  research  effort. 

2.1  Description  of  the  DTA 

In  its  present  state,  the  Dynamic  Test  Article  (Fig.  2. 1 )  is  composed  of  four  substructures: 
the  ring  truss,  two  solar  arrays,  and  an  equipment  tripod.  In  order  to  control  the  structure  using 
active  measures,  eight  mass  actuators  are  also  mounted  on  the  structure.  In  addition,  three 
mass  simulators  are  attached  to  the  structure  beneath  the  ring  truss.  These  mass  simulators 
replaced  the  original  DTA  box  truss,  which  could  not  be  relocated  to  AFIT  (3:15). 

2.1.1  Ring  Truss.  The  first  substructure  to  be  designed,  the  ring  truss  is  the 
structural  backbone  of  the  DTA.  The  circular  truss,  comprised  of  aluminum  tubes  bonded 
between  aluminum  joint  blocks,  is  designed  to  withstand  the  gravitational  load  of  the  testbed 
(7:41).  Six  reaction  mass  actuators  and  three  mass  simulators  are  mounted  on  the  ring  truss, 
simulating  the  presence  of  on-board  systems.  The  actuators  and  mass  simulators  were  mounted 
to  preserve  the  axis  of  symmetry  of  the  structure  (6:7). 

2.1 .2  Tripod.  The  DTA  tripod  consists  of  a  secondary  steel  plate,  three  rectangular 
aluminum  legs,  and  mounting  hardware  to  connect  the  tripod  to  the  ring  truss.  Designed  to 
have  low  frequency  flexible  modes,  the  tripod’s  design  was  based  on  a  Cassegrain  system  with 
the  tripod  supporting  the  secondary  reflecting  surface  (9:10).  Mounted  on  the  secondary  steel 
plate  are  two  reaction  mass  actuators.  The  actuators’  axes  are  perpendicular  to  each  other  and 
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Figure  2.1  DTA  Configuration  (6) 

parallel  to  the  ground.  The  legs  of  the  tripod  are  connected  to  the  ring  truss  by  aluminum 
mounting  plates. 

2.1 .3  Solar  Arrays.  Attached  to  the  ring  truss,  the  solar  arrays,  designed  to  simulate 
the  behavior  of  actual  solar  arrays,  are  sized  to  have  the  ability  to  provide  the  necessary  on- 
orbit  power  required  by  the  Cassegrain  optical  system.  Each  solar  array  is  comprised  of  a 
gridwork  of  thin  aluminum  strips  attached  to  a  15  ft.  aluminum  mast.  The  solar  arrays  were 
designed  to  simulate  the  low  frequency  behavior  of  actual  solar  arrays  as  those  modeled  for 
the  RSA.  The  solar  arrays  are  mounted  on  the  ring  truss  using  an  aluminum  block  assembly 
(9:61). 

2.1.4  Suspension  System.  The  DTA  is  supported  by  three  CS A  Engineering  Zero 
Spring  Rate  Mechanism  (ZSRM)  supports.  Figure  2.2  shows  the  relationship  between  the 
suspension  system  and  the  DTA  testbed.  The  ZSRMs,  which  are  attached  to  the  DTA  support 
frame,  are  able  to  support  the  testbed  by  supplying  air  pressure  beneath  a  piston  to  which 
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ZERO  SPRING 


Figure  2.2  DTA  and  Support  Structure  (6) 


the  support  cables  are  attached.  The  pistons  float  on  air  bearings  resulting  in  a  soft,  nearly 
ftictionless  suspension  (4:30).  Air  tanks  attached  to  the  support  frame  provide  air  for  each 
piston  (3:21).  The  support  cables  are  attached  to  the  ring  truss  in  three  locations  positioned  in 
order  to  maintain  the  axis  of  symmetry  of  the  testbed.  Thus,  the  testbed  is  isolated  from  any 
external  disturbances. 

2.2  Active  Vibration  Control  System 

An  active  vibration  suppression  system  is  incorporated  into  the  design  of  the  DTA.  The 
suppression  system  consists  of  eight  actuator/sensor  units  which  are  used  to  control  flexible 
body  modes.  Each  actuator/sensor  unit,  shown  in  Fig.  2.3,  is  comprised  of  a  reaction  mass 
actuator  (RMA)  and  a  Sundstrand  QA-1400  accelerometer.  Each  unit  produces  two  sensor 
signals:  (1)  inertial  acceleration  of  the  actuator  frame  at  the  mounting  location  and  (2)  the 
relative  velocity  between  the  reaction  mass  and  the  actuator  frame.  The  positioning  of  the 
actuator/sensor  units  is  shown  in  Fig.  2.1. 
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Figure  2.3  Actuator/Sensor  Unit 


22.1  Reaction  Mass  Actuator.  The  RMA  is  comprised  of  a  motor  coil,  a  linear 
velocity  transducer  (LVT),  a  reaction  (or  proof)  mass,  extension  springs,  and  a  housing.  For 
the  six  actuators  mounted  on  the  ring  truss,  two  extension  springs  are  mounted  to  both  the 
reaction  mass  and  the  assembly  housing  to  provide  for  gravity  off-load  for  the  mass.  The 
actuator  has  a  natural  frequency  of  1 .5-Hz  and  a  stroke  length  of  approximately  ±  1  in.  (9:90). 
The  two  remaining  actuators,  mounted  horizontally  on  the  secondary  steel  plate,  require  four 
extension  springs.  The  RMA,  which  is  attached  to  the  DTA  via  the  housing,  measures  the 
relative  velocity  between  the  motor  mass  and  actuator  housing.  Power  is  supplied  to  the 
actuators  via  the  PACOSS  Motor  Control  Units  (3:34-38). 

2.2.2  Accelerometer.  A  Sundstrand  QA-1400  accelerometer  is  attached  to  each 
actuator  housing  to  measure  the  inertial  acceleration  of  the  test  article  at  that  location.  If 
control  algorithms  were  to  be  implemented,  these  inertial  acceleration  measurements  would 
be  integrated  to  obtain  the  inertial  velocity  (3:34).  The  QA-1400  has  a  bandwidth  of  at  least 
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Figure  2.4  Motor  Control  Unit:  Front  Panel  (6) 

100  Hz;  therefore,  the  accelerometer  measurement  may  be  treated  as  a  perfect  measurement 
for  system  identification  purposes.  Signal  processing  (i.e.,  integration  of  the  measurement),  if 
desired,  can  be  supplied  by  the  PACOSS  Motor  Control  Units  (9:90). 

2.2.3  Motor  Control  Unit.  The  Motor  Control  Unit,  shown  in  Fig.  2.4,  provides  all 
of  the  signal  processing  and  power  for  an  actuator/sensor  unit.  There  are  a  total  of  eight  control 
units — one  for  each  actuator/sensor.  Each  unit  is  comprised  of  a  linear  motor  power  amplifier, 
an  accelerometer  pre-amp,  an  integrator,  and  an  LVT  amplifier,  and  adjustable  potentiometers. 
Adjustable  gain  controls  and  measurement  output  ports  are  included  on  the  front  panel  of 
the  instrument  cabinet  (9:90,93).  The  velocity  output  port  on  each  unit  was  connected  to  an 
analyzer  to  generate  frequency  and  time  response  data. 

2.3  Optima! 3 

Built  by  Systolic  Systems  Inc.,  the  Optima/3  Real-Time  Data  Acquisition  and  Control 
Processor  is  a  high-speed  parallel  processor  that  is  designed  for  data  acquisition,  estimation 
and  control  applications.  Comprised  of  a  Sun  3  master  processor,  eight  Systolic  parallel 
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processors,  and  a  data  acquisition  system,  the  Optima  has  16  input  channels  and  16  output 
channels.  The  Optima’s  main  purpose  in  this  research  effort  was  friction  compensation.  The 
sampling  rate,  signal  range,  and  input  gain  for  each  channel  were  programmed  specifically  for 
use  in  identifying  a  large,  flexible  space  structure  (10:2-1).  The  Optima/3  has  valid  sampling 
rates  ranging  from  2  to  20,000  Hz  with  a  2  Hz  resolution.  This  corresponds  to  sampling 
periods  of  0.05  to  500  msec  (10:3-4). 

Friction  compensation  was  accomplished  using  software  written  in  C  on  the  Sun  3/50 
host.  The  LABWARE^^  download  utility  was  able  to  download  the  compiled  C  code  to  the 
Optima/3.  The  Optima  then  executes  the  friction  compensation  program  (10:3-1). 

2.4  Additional  Hardware  and  Software 

A  Tektronix  2642A  Personal  Fourier  Analyzer  was  used  to  collect  the  data  and  to  create 
the  fi-equency  response  functions.  The  2642A  analyzer  is  a  high  performance  analog  data 
acquisition  and  signal  processing  peripheral  for  use  with  personal  computers.  The  analyzer 
processes  data  at  bandwidths  up  to  200  kHz,  with  a  sample  rate  2.56  times  the  bandwidth. 
The  analyzer,  however,  can  only  acquire  and  process  four  channels  of  data  (12:1-1).  The 
requirements  for  this  research  effort  was  a  minimum  of  8  input  channels.  Therefore,  a 
Tektronix  SI  5010  Programmable  Scanner  was  used  to  increase  the  number  of  available  input 
channels  from  four  to  sixteen. 

There  were  two  software  packages  that  were  used  to  collect  and  process  the  data.  To 
collect  the  data,  the  FSCAN  utility  (11)  and  the  Instrument  Program  (12)  from  Tektronix  were 
used  in  conjunction  with  the  Tektronix  2642A  Fourier  Analyzer.  The  FSCAN  utility,  running 
on  a  personal  computer,  interfaces  with  the  SI  5010  scanner  and  the  Instrument  Program  (IP)  to 
allow  data  to  be  taken  in  a  shorter  amount  of  time.  The  other  software  package,  MATLAB^^ 
from  MathWorks,  Inc.,  was  used  in  the  analysis  of  the  data  (13).  The  ERA  algorithm  was 
coded  in  MATLAB^^  and  was  used  to  generate  the  state-space  model  for  the  PACOSS  DTA. 
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III.  Theory 


The  material  presented  in  this  chapter  represents  the  theoretical  basis  for  the  develop¬ 
ment  of  an  accurate  state-space  model  which  can  be  used  to  design  a  vibration  suppression 
system.  First,  the  method  for  determining  the  amount  of  Coulomb  damping  present  in  the 
reaction  mass  actuators  will  be  developed.  Knowing  the  amount  of  Coulomb  damping  allows 
for  the  removal  of  this  nonlinearity  from  the  frequency  response  functions.  The  next  section 
explains  the  Eigensystem  Realization  Algorithm  which  will  be  used  to  generate  the  state-space 
model  from  the  frequency  response  functions. 

3.1  Coulomb  Damping 

The  reaction  mass  actuator,  acting  as  the  control  mechanism  for  the  testbed,  is  the  only 
moving  pan  on  the  DTA.  Like  every  other  mechanical  system,  friction  is  also  present  in  this 
mechanism.  In  this  case,  Coulomb  (dry  friction)  damping  is  the  major  friction  force  at  work. 
Coulomb  damping  cannot  be  physically  removed  from  the  system.  However,  by  employing  a 
nonlinear  compensator,  the  effects  of  Coulomb  damping  can  be  removed  from  the  frequency 
response  functions. 

The  RMA,  which  consists  of  two  extension  springs,  a  reaction  mass,  and  an  LVT,  is 
modeled  as  a  single-degree-of-freedom  harmonic  oscillator.  Figure  3.1  shows  a  schematic 
of  the  RMA.  The  Coulomb  damping  force,  Fc,  always  opposes  the  motion  of  the  mass. 
Therefore,  if  the  velocity  of  the  mass  is  negative,  then  the  friction  is  positive  and  vice  versa. 
The  equation  of  motion  for  the  RMA  is 

mx  -t-  Fc  sgn(x)  +  kx  =  0  (3.1) 


where 

sgn(i)  =  i| 

accounts  for  the  sign  change.  The  initial  conditions  are  a:(0)  =  xo  and  i;(0)  =  0. 
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Figure  3.1  Functional  Diagram  of  Reaction  Mass  Actuator 


The  piecewise  linear  solution  to  Eq.  (3.1)  is 

x{t)  =  A  cos  Unt  +  BsinuJnt  —  Fc  sgn{x)/k  (3.2) 

where  —  k/m  and  A  and  B  are  constants.  Equation  3.2  is  only  valid  between  velocity 
sign  changes.  Applying  the  initial  conditions,  the  solution  becomes 

x{t)  =  (  Fc  sgn{x)lk  +  xo)  cosa;„f  -  Fc  sgn(i;)  jk  (3.3) 

The  output  of  the  LVT,  which  is  the  response  that  will  be  used  to  determine  the  friction 
force,  Fc,  is  the  velocity  of  the  proof  mass  relative  to  the  actuator  housing.  Therefore,  the 
velocity  of  the  proof  mass  (x)  is  most  important.  Taking  the  derivative  of  x{t),  the  velocity 
function  is 

^(0  =  -{Fcsgn{x)/k  +  xo)oJn  sinuj  (3.4) 
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Figure  3.2  Free  Decay  Response  Function 

Shown  in  Fig.  3.2  is  a  representative  free  decay  response  for  both  displacement  and 
velocity  of  the  actuator.  To  solve  for  the  Coulomb  damping  force  present  in  the  mechanism, 
piecewise  techniques  will  be  used.  The  two  regions  highlighted  in  Fig.  3.2  will  be  considered. 
The  first  region  is  where  x  is  negative;  the  second,  where  x  is  positive. 

3.1.1  Case  1:  x  <0.  Throughout  Region  I  (with  the  exception  of  the  start  and  end 
times),  the  velocity  is  negative.  Substituting  this  condition  into  Eqs.  (3.3)  and  (3.4)  yields 

xi{t)  =  (a:o  -  Fc  /k)  cosu;„f  +  Fc  /k  (3.5a) 

=  -{^0-  Fc  /k)u}n  s'mujnt  (3.5b) 

Now,  xi{ti)  and  must  be  calculated,  where  the  subscript  I  denotes  the  responses  in 
Region  I,  and  ti  is  the  time  in  Region  I  where  i  =  0.  From  Eq.  (3.5b),  for  x  to  equal  0  then 

sin  u>nt  =  0 
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Thus,  =  TT  /  Un.  Substituting  t  =  Uin  Eqs.  (3.5a)  and  (3.5b)  yields 


xi{ti)  =  2  Fc  I k  -  xo  (3.6a) 

ii(^i)  =  -(a^o  -  f'c  /k)u!n  sin  ujnh  =  0  (3.6b) 

Referencing  Fig.  3.2,  ti  is  the  end  time  for  Region  I  and  the  start  time  for  Region  II;  thus, 
a;i(fi)  and  xi(^i)  become  the  initial  conditions  for  Region  II. 

3.1.2  Case  II:  i  >  0.  In  Region  II,  with  the  exception  of  the  start  and  end  times 
once  again,  the  velocity  is  positive.  From  above,  the  initial  conditions  for  this  region  are 

a:ii(fi)  =  2Fc  jk  -  xo 

and 

=  0 

Using  trigonometric  formulas  and  the  condition  i  >  0,  the  solutions  in  this  region  can  be 
rewritten  as 

a;ii(f)  =  Asinun{t  —  fi)  +  B cos u>n{t  —  ti)  —  Fc  Ik  (3.7a) 

i:il(f)  =  A  U>n  COSUJn{t  —  t])  —  BuJn  sinUn{t  —  ti)  (3.7b) 

Applying  the  initial  conditions  for  Region  H,  Eq.  (3.7)  becomes,  for  ti  <  t  <t2, 

xii{t)  =  {3  Fc  Ik  -  xo)  cos  u:n{t  -  fi)  -  Fc  Ik  (3.8a) 

in{t)  =  -{3  Fc  Ik  -  xo)uin  sina;„(f  -  fi)  (3.8b) 
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Next,  t2  must  be  determined.  ^2  is  the  next  zero  crossing  of  the  velocity  function  in  the  region. 
Again  referencing  Fig.  3.2,  ^2  is  found  to  be 

TT 

^2  =  H - 


Evaluating  Eq.  (3.8)  at  t  =  t2  yields 


a;ii(^2)  =  xq-  A  Fc  Ik  (3.9a) 

xn{t2)  =  -(3Fc  /k  -  xo)ujn  sinu;„(— )  =  0  (3.9b) 

3.1 .3  Determining  the  Coulomb  Damping  Coefficient.  Using  the  peak  amplitudes 
of  Eqs.  (3.6b)  and  (3.9b),  the  decrease  in  amplitude  from  t  =  ti  to  t  =  t2  is  calculated.  The 
decrease  over  every  half  cycle  is  found  to  be 


(xo  -  Fc  Ik)  Un  -  {xo  -SFc  lk)LOn  =  {2  Fc  I k)  Un  (3.10) 


Thus,  the  amplitude  decays  by  (4  Fc  lk)ujn  every  cycle. 

Recalling  that  ojn  =  2wfn  where  /„  is  the  natural  frequency  in  Hz,  the  decrease  in 
amplitude  per  cycle.  Ay,  can  be  written  as 


A  ^Fctt 

A!,  =  — /. 


Rearranging,  the  Coulomb  damping  force  can  be  determined  by 


Fn  = 


Ay  k 

StT  fn 


(3.11) 


In  all  probability,  more  than  one  cycle  will  be  used  to  determine  the  coefficient;  thus,  the  more 
useful  formula  is 


Ay  k 

8  TT  /„  (no.  of  cycles) 


(3.12) 
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where  Ay  is  redefined  as  the  decrease  in  amplitude  over  the  number  of  cycles  to  be  used  in 
determining  the  Coulomb  damping  coefficient. 

Using  Eq.  (3.12),  the  amount  of  Coulomb  damping  in  the  actuator  can  be  calculated.  This 
is  of  importance  since  the  friction  force  can  be  mathematically  removed  from  the  frequency 
response  functions.  In  removing  this  nonlinearity  from  the  frequency  response  functions,  it  is 
hoped  that  a  more  accurate  state-space  representation  of  the  testbed  will  be  generated. 

3.2  System  Identification 

The  main  goal  of  this  thesis  is  to  obtain  an  accurate  mathematical  representation  of  the 
PACOSS  DTA  from  frequency  response  functions.  The  basic  problem  posed  is  one  of  system 
identification:  Given  measurement  functions  y{k),  construct  a  state-space  realization  such 
that  the  functions  y  are  reproduced  by  the  state-variable  equations  (5:621). 

ERA  is  based  upon  a  singular  value  decomposition  of  the  block  Hankel  matrix,  which  is 
composed  of  impulse  response  data.  This  modal  synthesis  technique  is  capable  of  identifying 
the  system  if  there  is  little  or  no  noise  found  in  the  measurements  (2:1564). 

Consider  the  discrete-time  linear  system: 

x{k  +  1)  =  Ax{k)  +  Bu{k)  (3.13a) 

y(^k)  =  Cx{k)  +  Du{k)  (3.13b) 

where  x  is  an  n-dimensional  state  vector,  u  an  p-dimensional  input  vector,  and  y  a  y- 
dimensional  output  vector.  The  A,  B,C,  and  D  matrices  have  dimensions  (nxn),(nxp),(y 
X  n),  and  (qxp),  respectively.  The  integer  k  is  the  sample  indicator.  The  A  matrix  represents 
the  dynamics  of  the  system.  For  flexible  structures,  it  is  a  representation  of  the  mass,  stiffness, 
and  damping  properties  (5:621).  For  the  system  in  Eq.  (3.13)  with  unit  impulse  response,  the 
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time-domain  description  is  given  by  the  function  known  as  the  Markov  Parameters  ^ : 

Y{k)  =  (3.14) 

The  ERA  technique  first  forms  a  block  Hankel  matrix.  As  mentioned  before,  this  matrix  is 
composed  of  the  sampled  unit  impulse  response: 

■  Y{k)  •••  Y{k  +  ms-i)  - 
H{k-l)=  !  ••.  :  (3.15) 

_Y[k-\-lr-i)  •••  K(A;  4- /r-i  +  m^-i) _ 

where  r  and  5  are  arbitrary  integers  satisfying  the  inequalities  rq  >  n  and  sp  >  n,  and 
li  {i  =  1, 2, . . . ,  r  —  1)  and  rrij  {j  =  1, 2, . . . ,  s  —  1)  are  arbitrary  integers.  The  singular 
value  decomposition  of  the  r  x  5  block  Hankel  matrix  is  expressed  as: 

/f(0)  =  PDQ'^  (3.16) 

Using  the  decomposition  of  the  Hankel  matrix,  the  state-space  realization,  which  is  discrete¬ 
time  and  reduced-order,  is  computed: 

A  =  (3.17a) 

B  =  D]!^QIE,  (3.17b) 

C  =  ElPnD]!^  (3.17c) 

D  =  r(0)  (3.17d) 

The  matrices  Pn  and  Qn  are  composed  of  the  first  n  columns  of  P  and  Q  from  the  decompo¬ 
sition.  Dn  is  a  diagonal  matrix  of  n  non-zero  singular  values.  E^  is  defined  as  [Ip,  0]  and  Ej 

^In  literature  the  discrete-impulse  response  function  sequence  is  commonly  referred  to  as  the  Markov 
Parameters. 
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is  defined  as  [/„  0],  where  Ip  and  Iq  are  identity  matrices  of  the  order  p  and  q,  respectively, 
and  0  is  the  zero  matrix  (2:1565). 

For  a  more  in-depth  look  at  the  Eigensystem  Realization  Algorithm,  see  the  original 
derivation  by  Juang  and  Pappa  (5:620-627). 
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IV.  Experimental  Procedure 


The  following  sections  detail  the  procedures  used  to  obtain  and  process  the  data  in 
this  research  effort.  Procedures  used  to  identify  the  Coulomb  damping  coefficients  and 
to  implement  the  compensation  routine  will  be  discussed.  The  procedures  that  were  used 
to  acquire  the  frequency  response  data,  compensated  and  non-compensated,  also  will  be 
presented.  The  method  used  to  validate  the  ERA  routines  will  also  be  discussed. 

4.1  DTA  and  Actuator  Preparation 

In  1992,  Capt  Scott  George  completed  a  finite  element  analysis  of  the  DTA.  To  obtain 
the  finite  element  model  of  the  DTA,  151  accelerometers  and  the  associated  instrumentation 
wiring  were  connected  to  the  DTA  testbed  (4:73).  The  accelerometers  and  wires  remained 
attached  to  the  structure  throughout  the  next  two  research  efforts.  Due  to  the  weight  and 
stiffness  of  the  wires,  damping  and  stiffness  were  being  added  to  the  frequency  response  data. 
All  of  the  accelerometers  and  the  instrumentation  wiring  were  removed  from  the  structure, 
since  they  were  no  longer  needed.  Thus,  the  damping  and  stiffness  that  resulted  from  the 
wires  was  removed  from  the  frequency  response  functions. 

Maintenance  was  also  performed  on  the  eight  mass  actuators.  Due  to  the  fact  that  the 
masses  had  been  suspended  by  the  same  springs  since  at  least  late  199T;  the  extension  springs 
on  the  eight  RMAs  were  replaced.  Upon  replacement,  the  net  spring  constant  was  calculated 
via  bench  testing.  The  net  spring  constant,  k,  was  found  to  be  equal  to  1.0  Ib/in.  This 
verified  the  results  obtained  by  Martin  Marietta  (9:90).  The  actuator’s  shaft,  linear  velocity 
transducer,  and  linear  bearings  were  cleaned  and  lubricated  with  a  silicone-based  lubricant  to 
eliminate  friction  in  the  system.  This  resulted  in  improved  results  when  free  decay  tests  were 
accomplished  to  determine  the  Coulomb  damping  coefficients. 
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Y1  Y2  Y3  Y4 


Figure  4.1  Truth  Model  Diagram 
4.2  Validation  of  ERA  Routines 

The  ERA  routines  which  performed  the  system  identification  were  written  in  MATLAB^^ 
by  Capt.  Richard  Cobb,  a  Ph.D.  candidate  at  the  Air  Force  Institute  of  Technology  (1).  Several 
of  the  routines  were  modified  for  use  in  this  thesis.  Before  using  these  routines  to  identify  the 
PACOSS  DTA,  a  smaller  scale  identification  was  used  to  validate  the  usefulness  of  the  code. 

The  spring-mass-damper  system  pictured  in  Fig.  4. 1  was  used  as  the  truth  model  for  the 
validation  of  the  ERA  routines.  This  system  was  modeled  after  the  system  simulator  box  built 
by  Tektronix  for  use  with  the  Instrument  Program  tutorial  (12:1-9).  The  system  is  composed 
of  four  equal  masses  (1  x  10“^  slugs),  interconnected  with  four  springs  and  five  dampers. 
The  four  springs  have  the  same  spring  constant  (1  x  10®  Ib/in)  and  the  dampers  have  the 
same  damping  coefficient  (2  lb- sec/in).  There  are  fom  force  inputs  (FI  through  F4)  and  four 
displacement  outputs  (Y1  through  Y4)  in  the  system. 

Three  inputs  (FI,  F2,  and  F3)  and  three  outputs  (Yl,  Y2,  and  Y3)  were  used  to  obtain 
a  3  x  3  MIMO  system.  Two  cases  were  used:  (1)  the  truth  model  without  noise  and  (2)  the 
truth  model  with  a  15%  measurement  noise.  Noise  was  added  to  the  system  to  ascertain  the 
accuracy  of  the  algorithm  without  perfect  measurements.  The  noise  was  added  directly  to 
the  FRFs.  This  is  not  the  case,  however,  with  the  PACOSS  DTA,  where  noise  is  continually 
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introduced  into  the  system.  The  rms  value  of  the  added  noise  was  obtained  by  taking  15%  of 
the  rms  value  of  the  frequency  response  function.  This  noise  was  then  added  to  the  FRFs. 

Frequency  response  functions  for  the  truth  model  system  were  generated  using  MATLAB^^ . 
The  equations  of  motion  and  a  state-space  representation  of  the  system  were  derived  from 
Fig.  4.1.  The  state-space  representation  was  discretized  and  the  frequency  response  functions 
were  calculated  using  mkfrf.m.  The  15%  noise  was  added  using  addnoise.m  (1). 

The  ERA  routines  were  then  employed  to  find  a  state-space  realization  for  the  truth 
model,  with  and  without  noise.  The  parameter  settings  for  the  system  identification  can  be 
found  in  Appendix  J.  The  undamped  natural  frequency  and  damping  ratio  of  the  poles  of  the 
ERA  fit  for  both  systems  will  be  compared  with  the  poles  of  the  truth  model.  This  comparison 
will  be  used  as  the  metric  to  determine  the  accuracy  of  the  ERA  fits. 

4 .3  Coulomb  Damping  Identification 

As  mentioned  in  Chapter  HI,  free  decay  was  used  to  determine  the  amount  of  Coulomb 
damping  in  the  eight  reaction  mass  actuators,  which  provide  the  disturbances  to  generate  the 
frequency  response  functions  (14:125).  In  free  decay,  the  mass  is  displaced  a  fixed  distance 
and  released.  The  spring-mass  system,  shown  in  Fig.  3.1,  will  oscillate  until  the  total  energy 
left  in  the  system  can  no  longer  overcome  gravity  and  the  friction  forces  present.  The  data 
obtained  from  this  method  consists  of  a  voltage  proportional  to  the  velocity  of  the  proof  mass 
relative  to  the  actuator  housing.  The  voltages,  obtained  from  the  LVT  OUT  port  on  the  motor 
control  unit,  were  stored  by  the  IP  program  and  converted  to  ASCII  format. 

A  MATLAB^^  routine,  which  is  referred  to  as  a  m-file,  was  written  to  analyze  the  time 
response  data.  The  m-file  (shown  in  Appendix  L)  plots  the  data  and  calculates  the  frequency, 
the  number  of  cycles,  and  the  change  in  amplitude  during  those  cycles.  The  routine  uses 
Eq.  (3.12)  and  the  free  decay  data  to  calculate  the  Coulomb  damping  coefficient.  The  net 
spring  constant,  k,  is  equal  to  1  Ib/in  (9:90). 

For  each  actuator,  four  sets  of  data  were  taken.  Coulomb  damping  coefficients  were 
calculated  from  each  set  of  data.  The  four  coefficients  obtained  for  each  actuator  were  then 
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averaged  to  find  an  approximation  of  the  Coulomb  damping  coefficient  for  that  actuator.  To 
determine  the  accuracy  of  the  Coulomb  damping  coefficients,  a  Coulomb  damping  compen¬ 
sator,  which  was  programmed  in  C  code  and  downloaded  to  the  Optima/3,  was  used.  If  the 
coefficients  used  are  accurate,  the  compensation  force  applied  should  not  induced  any  oscil¬ 
latory  motion  since  the  force  applied  should  be  the  exact  amount  to  neutralize  the  Coulomb 
damping  force.  In  a  test  to  determine  the  accuracy  of  the  calculated  coefficients,  the  com¬ 
pensation  force  was  the  only  force  applied.  The  compensation  force  caused  each  mass  to 
oscillate;  thus,  thus,  the  calculated  coefficients  were  too  large  and,  if  used,  would  result  in  an 
unstable  system.  A  method  was  then  devised  to  adjust  the  Coulomb  damping  coefficient  for 
each  actuator. 

Using  the  compensator,  Fourier  analyzer  and  IP  program,  free  decay  was  employed  to 
adjust  the  damping  coefficients.  When  the  correct  coefficient  was  used,  there  was  no  decrease 
in  the  amplitude  of  the  time  response  function.  Viewing  the  free  decay  response  via  the 
IP  program,  the  coefficient  could  be  adjusted  either  up  or  down  to  provide  the  desired  time 
response.  The  C  code  used  to  adjust  these  coefficients  can  be  found  in  Appendbc  K.  The 
adjusted  coefficients  where  then  programmed  into  the  Coulomb  damping  compensator  and 
used  to  remove  the  effects  of  Coulomb  damping  from  the  frequency  response  functions. 

4.4  Coulomb  Damping  Compensation 

Coulomb  damping  compensation  was  accomplished  by  using  the  Optima/3.  A  nonlinear 
compensator  was  programmed  in  C  code  to  mathematically  remove  the  effects  of  Coulomb 
damping  from  the  frequency  response  functions. 

The  couldamp.c  program  can  be  found  in  Appendix  K.  The  compensator  sampled  eight 
relative  velocities  for  the  eight  mass  actuators  200  times  per  second.  The  velocities,  which 
were  output  by  the  LVTs  and  read  in  by  the  Optima  as  voltages,  were  converted  to  inches 
per  second.  The  sign,  as  well  as  the  magnitude,  of  the  velocity  was  checked.  This  is  due 
to  the  signum  function  in  Eq.  (3.4).  However,  a  “dead  zone”  was  created  at  ±0.2  inches  per 
second.  The  dead  zone  was  created  because  the  corresponding  voltages  for  velocities  below 
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0.2  inches/second  are  the  same  order  of  magnitude  as  noise  found  in  the  system.  Without  this 
dead  zone,  the  compensator  would  apply  Coulomb  damping  compensation  when  not  needed. 
It  should  also  be  mentioned  that  85%  of  the  adjusted  Coulomb  damping  coefficient  for  each 
actuator  was  used.  This  was  done  so  that  the  actuators  were  not  overcompensated,  which 
would  yield  an  unstable  system. 

If  the  mass’  relative  velocity  was  greater  than  0.2  inches  per  second,  then  a  positive 
friction  compensation  was  applied  by  sending  a  voltage  proportional  to  the  friction  coefficient. 
For  a  negative  velocity,  a  negative  compensation  voltage  is  applied.  The  Optima  added  the 
compensation  voltage  to  the  random  noise  voltage  that  was  output  by  the  Tektronix  2642A 
analyzer.  This  sum  was  then  output  to  the  motor  control  unit  of  the  excitation  actuator. 

4 .5  Data  Acquisition 

4.5.1  Frequency  Response  Functions.  The  state-space  representation  of  the  DTA 
was  identified  based  only  on  experimental  measurements  of  the  open  loop  system.  Figure  4.2 
shows  the  equipment  configuration  used  to  obtain  the  frequency  response  functions  (FRFs). 


Figure  4.2  Equipment  Configuration  for  Open  Loop  System 


The  FRFs  were  generated  by  the  2642A  Fourier  analyzer.  The  analyzer,  outputting  a 
600  mV  random  noise  with  a  20  Hz  bandwidth,  was  used  to  excite  an  actuator.  The  analyzer 
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Figure  4.3  Equipment  Configuration  with  Active  Friction  Compensation 

and  IP  program,  in  conjunction  with  the  SI  5010  scanner,  obtained  the  inertial  acceleration 
readings  from  each  accelerometer.  The  analyzer  used  this  data  to  calculate  the  corresponding 
FRFs.  This  was  done  seven  more  times  to  obtain  a  complete  set  of  64  FRFs.  The  accelerometer 
readings  were  sampled  2048  times  per  frame  with  a  20  Hz  base  bandwidth  using  a  Hanning 
windowing  type  with  no  overlap.  Hanning  windowing,  which  is  recommended  for  structural 
vibration  data,  suppresses  leakage  of  the  data.  These  settings  give  a  sampling  rate  of  51.2 
samples  per  second  and  a  frequency  resolution  of  0.025  Hz  (12:A-16). 

4.5.2  Frequency  Response  Functions  with  Friction  Compensation.  As  mentioned 
earlier,  the  Coulomb  friction  force  in  the  actuators  introduces  nonlinearities  into  the  system. 
Therefore,  active  friction  compensation  was  employed  to  remove  the  effect  of  this  friction 
force  from  the  FRFs.  Figure  4.3  illustrates  the  equipment  configuration  used  to  obtain  this  set 
of  frequency  response  data. 

The  IP  program,  scanner,  and  analyzer  settings  for  this  setup  were  virtually  identical. 
The  only  difference  was  that  a  300  mV  random  noise  was  output  from  the  analyzer  to  act 
as  a  disturbance.  The  output  voltage  level  was  decreased  to  ensure  that,  with  active  friction 
compensation,  the  actuator  did  not  saturate  and  hit  the  stops. 
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V;  Results 


The  following  sections  describe  the  results  of  the  procedures  previously  discussed  in 
Chapter  IV.  The  discussion  includes  the  calculated  and  adjusted  Coulomb  damping  coefficients, 
comparison  of  truth  model  FRFs  and  the  corresponding  ERA  fit,  results  of  the  ERA  algorithm 
when  applied  to  the  FRFs  obtained  from  the  PACOSS  DTA,  and  results  of  a  comparison 
between  FRFs  with  and  without  active  friction  compensation  as  well  as  FRFs  obtained  by 
Lt.  Matheson  in  1992.  It  should  be  noted  that  the  measured  data  was  acceleration  frequency 
response  functions.  One  integration  was  performed  on  the  data  to  obtain  velocity  frequency 
response  functions.  The  results  shown  in  this  chapter  are  the  ERA  fits  to  the  velocity  FRFs. 

Throughout  this  chapter,  the  terms  collocated  and  non-collocated  are  used.  Collocated 
is  used  to  describe  frequency  response  functions  when  the  measurement  is  obtained  from  the 
same  point  where  the  excitation  force  is  applied.  For  example,  a  frequency  response  function 
for  Measurement  1  Excitation  Actuator  1  is  a  collocated  FRF  since  the  measurement  was 
taken  at  Actuator  1  and  the  force  was  also  applied  there.  Non-collocated  is  used  to  describe 
FRFs  when  the  measurement  is  obtained  at  a  different  location  than  where  the  excitation  force 
was  applied.  An  example  of  a  non-collocated  FRF  is  Measurement  1  Excitation  Actuator  2. 

5.1  Coulomb  Damping  Coefficients 

Table  5 . 1  shows  the  Coulomb  damping  coefficients  that  were  derived  experimentally  and 
the  adjusted  coefficients,  which  were  used  in  the  compensation  program.  During  the  course 
of  the  research,  it  was  found  that  these  coefficients  change  depending  on  the  temperature  and 
humidity. 

As  shown  in  Table  5.1,  there  is  a  big  discrepancy  in  the  coefficients  for  Actuators  9 
and  10.  This  is  explained  by  the  fact  that  these  two  actuators  are  mounted  horizontally.  The 
linear  bearings  used  in  the  eight  actuators  are  not  designed  for  horizontal  use.  Thus,  when 
free  decay  was  used,  only  one  cycle  occurred  before  the  motion  died  out  resulting  in  a  rather 
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Table  5 . 1  Coulomb  Damping  Coefficients 


Actuator 

No. 

Calculated 

Coefficient 

(Ibf) 

Adjusted 

Coefficient 

(Ibf) 

1 

0.0591 

0.0325 

2 

0.0565 

0.0500 

3 

0.0549 

0.0250 

4 

0.0353 

0.0100 

5 

0.0574 

0.0150 

6 

0.0518 

0.0500 

9 

0.2824 

0.0325 

10 

0.3188 

0.0175 

large  coefficient.  For  this  reason,  the  actuators  were  adjusted.  The  more  accurate,  adjusted 
coefficients  were  used  in  the  Coulomb  damping  compensator. 


5.2  ERA  Fit  of  the  Truth  Model 

As  mentioned  in  Chapter  IV,  a  smaller  MIMO  system,  the  truth  model  shown  in  Fig.  4. 1 , 
was  identified  first  to  validate  the  ERA  routines  that  were  used  in  the  identification  of  the 
PACOSS  DTA.  A  3  X  3  MIMO  system  utilizing /orce  inputs  FI  through  F3  and  displacement 
outputs  Y1  through  Y3,  with  and  without  noise,  was  used  as  the  validation  model. 

5.2.1  Truth  Model  without  Noise.  Shown  in  Fig.  5.1  and  5.2  are  two  FRFs  and 
their  corresponding  ERA  fits.  Figure  5.1,  the  transfer  function  and  the  corresponding  fit  for 
displacement  Y1  and  force  input  F2  is  representative  of  the  remaining  non-collocated  output 
measurement  and  force  input  pairs.  Figure  5.2  shows  the  collocated  transfer  function  and 
its  corresponding  fit  for  displacement  Y2  and  force  input  F2.  The  ERA  fit  for  the  collocated 
transfer  function  shown  in  Fig.  5.2  is  also  representative  of  the  remaining  collocated  transfer 
functions.  As  shown  in  Table  5.2,  the  algorithm,  in  the  case  with  no  noise,  is  able  to  accurately 
identify  the  system. 


5-2 


a  1 0OO  2000  3000  4000  5000  6000  7000  8000  9000  1 0OOO 

Frequency  (Hz) 


Figure  5.1  Displacement  Y1  Force  Input  F2  for  the  Truth  Model  with  No  Noise 


Figure  5.2  Displacement  Y2  Force  Input  F2  for  the  Truth  Model  with  No  Noise 
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Figure  5.3  Displacement  Y1  Force  Input  F2  for  the  Truth  Model  with  Noise 


Figure  5.4  Displacement  Y2  Force  Input  F2  for  the  Truth  Model  with  Noise 
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5.2.2  Truth  Model  with  Noise.  To  determine  how  well  the  Eigensystem  Realization 
Algorithm  could  handle  noise,  a  15%  noise  was  added  to  the  frequency  response  functions 
obtained  from  MATLABtm.  Figure  5.3  shows  the  transfer  function  and  the  corresponding 
ERA  fit  for  displacement  Y1  and  force  input  F2,  while  Fig.  5.4  shows  the  transfer  function 
for  displacement  Y2  and  input  F2  and  its  corresponding  fit.  The  ERA  fits  shown  for  these 
transfer  functions  are  characteristic  of  the  remaining  frequency  response  functions.  Even  in 
the  presence  of  noise,  the  ERA  algorithm  is  able  to  accurately  identify  the  truth  model.  As 
shown  in  Table  5.2,  the  ERA  algorithm  is  able  to  obtain  an  accurate  model  of  the  system  even 
with  15%  noise  included  in  the  frequency  response  functions. 


Table  5.2  Truth  Model  Results 


TRUTH  MODEL 

ERA  FIT 
with  No  Noise 

ERA  FIT 
with  Noise 

Mode 

Frequency 

c 

Frequency 

c 

Frequency 

c 

No. 

(Hz) 

(%) 

(Hz) 

(%) 

(Hz) 

(%) 

n 

1748.64 

5.02 

1748.64 

5.02 

1747.70 

4.98 

mm 

5031.80 

4.22 

5031.80 

4.22 

5031.66 

4.28 

H 

7709.75 

5.22 

7709.77 

5.22 

7708.31 

5.31 

H 

9458.39 

6.03 

9458.44 

6.03 

9457.83 

6.03 

53  ERAFitofthePACOSSDTA 

Based  upon  the  results  of  the  ERA  fit  on  the  truth  model,  it  was  determined  that  the 
Eigensystem  Realization  Algorithm  could  identify  the  poles  and  zeros  of  the  larger  system. 
The  state-space  representation  of  the  PACOSS  DTA  obtained  using  the  friction-compensated 
system  is  not  shown  due  to  the  size  of  the  system.  80  states  were  identified  under  14  Hz.  The 
resulting  A,  B,  C,  and  D  matrices  were  80  x  80, 80  x  8, 8  x  80,  and  8x8,  respectively.  The 
resulting  D  matrix  is  a  zero  matrix. 


The  results  presented  in  this  section  are  the  result  of  identifying  the  PACOSS  DTA 
with  active  friction  compensation.  All  transfer  functions  shown  are  for  inertial  acceleration 
outputs.  Figure  5.5  shows  the  results  for  a  non-collocated  sensor  and  excitation  actuator.  For 


this  particular  FRF,  the  measurement  was  taken  at  Actuator  1  with  Actuator  2  acting  as  the 
disturbance  input.  This  fit  is  representative  of  the  non-collocated  FRFs.  Figure  5.6  shows  the 
ERA  fit  for  a  collocated  sensor  and  actuator.  A  trend  was  noticed  in  the  identification  of  this 
structure:  the  ERA  fit  for  the  non-collocated  sensor/actuator  pairs  are  better  than  those  of  the 
collocated  pairs.  Figures  5.5  and  5.6  show  two  of  the  better  fits  that  were  obtained  from  the 
ERA  algorithm.  Not  all  of  the  ERA  fits  were  of  this  quality.  Figure  5.7  shows  the  frequency 
response  function  for  Measurement  1  Excitation  Actuator  6  and  its  corresponding  ERA  fit. 
This  FRF/ERA  fit  was  the  worst  in  the  set  of  64  transfer  functions;  however,  it  is  should  be 
noted  that  it  is  not  possible  to  fit  every  FRF  exactly  since  a  nonlinear  system  is  being  identified 
using  a  linear  scheme.  The  resulting  64  FRFs  and  their  corresponding  ERA  fits  are  shown  in 
Appendices  A-H. 


Figure  5.7  Measurement  1  Excitation  Actuator  6 
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5.4  Friction  Compensation:  Is  it  worth  it? 

One  of  the  goals  of  this  thesis  was  to  calculate  the  Coulomb  damping  coefficients  and 
apply  active  friction  compensation  to  the  PACOSS  DTA.  This  goal  was  achieved.  However, 
the  question  remains:  Was  friction  compensation  worth  the  trouble  to  calculate  the  coefficients 
and  implement  the  compensation  program?  The  answer  is  yes.  The  following  paragraphs 
give  a  more  in-depth  look  at  why  friction  compensation  should  be  used. 

First,  however,  the  measured  FRFs  of  both  the  compensated  and  non-compensated 
system  are  compared.  FRFs  from  each  system  are  shown  in  Figs.  5.8  and  5.9.  Figure  5.8 
compares  both  FRFs  for  a  non-collocated  sensor/actuator  parr  while  Fig.  5.9  compares  two 
collocated  FRFs. 


Figure  5.8  FRFs  for  Measurement  1  Actuator  2 

Referencing  both  of  these  figures,  one  can  see  that  the  only  difference  between  the  FRF 
of  the  compensated  system  and  the  FRF  of  the  non-compensated  system  is  the  actuator  mode 
at  1.5  Hz.  Of  course,  this  is  expected  since  the  friction  compensation  removes  the  nonlinear 
damping  from  that  mode.  What  is  surprising,  however,  is  that  the  actuator  mode  is  the  only 
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Figure  5.9  FRFs  for  Measurement  4  Actuator  4 

mode  that  changes.  Prior  to  this  comparison,  it  was  thought  that  many  of  the  low  frequency 
modes  would  be  altered  by  the  friction  compensation.  This  was  not  the  case. 

Due  to  the  fact  that  both  systems  resulted  in  almost  identical  FRFs,  it  was  assumed  that 
friction  compensation  was  not  worth  the  effort  and  that  the  ERA  fit  for  the  non-compensated 
frequency  response  functions  would  be  the  same  as  those  that  contained  the  friction  compen¬ 
sation.  Again,  this  was  not  the  case.  Looking  at  Figs.  5.10  and  5.1 1,  one  can  see  that  the  ERA 
fit  for  the  compensated  system  did  an  overall  better  job  identifying  the  poles  and  zeros  of  the 
DTA.  This  can  be  seen  in  the  identification  of  the  zeros  especially.  The  zeros  identified  in  the 
compensated  system  are  deeper  than  those  found  in  the  non-compensated  system.  Overall, 
the  ERA  fit  for  the  compensated  system  is  better  than  the  fit  for  the  non-compensated  system. 

A  reason  for  the  difference  in  the  ERA  fits  is  that  the  ERA  algorithm  produces  a  state- 
space,  or  linear,  model.  Thus,  nonlinearities  in  a  system  produce  larger  errors  in  the  fit. 
Removing  nonlinearities,  such  as  Coulomb  dtunping,  helps  to  produce  a  linear  system  to  be 
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Frequency  <Hz) 


Figure  5.10  ERA  Fit  for  Compensated  System:  Measurement  1  Actuator  1 


Figure  5.1 1  ERA  Fit  for  Non-Compensated  System:  Measurement  1  Actuator  1 
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used  for  system  identification.  Using  the  compensated  system,  then,  produces  a  “more  linear” 
model  which  is  more  accurately  identified  by  the  ERA  algorithm. 

5.5  Ejfects  of  Instrumentation  Wire  Removal 

In  1992,  Lt.  Chad  Matheson,  using  transfer  function  matching,  constructed  a  model  of 
the  PACOSS  DTA.  However,  the  test  article  at  that  time  was  outfitted  with  instrumentation 
wiring  for  151  accelerometers.  In  order  to  reduce  damping  and  the  number  of  nonlinearities 
in  the  system,  the  wiring  and  the  sensors  were  removed  so  that  an  accurate  model  could  be 
obtained.  After  obtaining  the  frequency  response  data,  it  was  compared  to  the  frequency 
response  data  obtained  by  Lt.  Chad  Matheson.  Figure  5.12  shows  the  frequency  response 
function  for  Measurement  2  Excitation  Actuator  1  obtained  from  the  compensated  system  and 
that  obtained  by  Lt.  Matheson  for  the  same  sensor/actuator  pair. 


Figure  5.12  Comparison  of  FRFs 

Up  to  9  Hz,  the  two  response  functions  match  rather  weU.  Starting  at  9  Hz,  the 
poles  and  zeros  of  the  two  systems  no  longer  match.  The  poles  and  zeros  not  only  occur  at 
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different  frequencies,  but  their  amplitudes  also  vary.  Figure  5.12  shows  that  the  removal  of 
the  instrumentation  wires  has  a  pronounced  effect  upon  the  system.  The  frequency  response 
functions,  both  compensated  and  non-compensated,  obtained  in  the  course  of  this  research  are 
reproducible  and  have  been  reproduced  several  times. 
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VI.  Recommendations  and  Conclusions 

This  thesis  had  two  primary  goals:  (1)  Obtain  an  accurate  state-space  model  of  the 
PACOSS  Dynamic  Test  Article  and  (2)  Identify  the  Coulomb  damping  in  each  actuator  and 
compensate  for  this  nonlinearity  in  the  frequency  response  functions.  The  two  primary  goals 
of  this  thesis  have  been  achieved.  The  state-space  representation  of  the  system’s  dynamics 
for  the  first  40  modes  is  accurate.  As  seen  in  Chapter  V,  the  Coulomb  damping  factor  in  each 
actuator  was  identified,  and  the  friction  compensation  was  successful,  as  demonstrated  by  the 
undamped  actuator  mode  in  the  compensated  system. 

As  demonstrated,  the  ERA  algorithm  worked  well  considering  the  dynamics  of  the 
structure  being  identified  is  nonlinear.  It  was  noticed  that  the  system  identified  by  ERA  did 
not  fit  the  collocated  sensor/actuator  pairs  as  well  as  the  non-collocated  pairs.  The  reason  for 
this  discrepancy  is  not  known  and  needs  to  be  investigated  further.  Overall,  the  ERA  algorithm 
did  a  good  job  of  identifying  the  poles  in  the  system,  which  is  important  for  stability  reasons. 

The  ERA  algorithm  performed  better  on  the  friction-compensated  system.  This  is 
expected  since  the  compensated  system  is  “more  linear”  than  the  non-compensated  system 
due  to  the  removal  of  the  damped  actuator  mode.  Thus,  the  friction-compensated  state-space 
model  should  be  used  for  control  design.  The  friction  compensation  can  be  taken  into  account 
in  the  design  of  the  controller.  This  makes  the  application  of  optimal  control  strategies  simpler 
since  the  structure’s  dynamics  can  be  considered  linear. 

As  depicted  in  Chapter  V,  the  removal  of  outside  influences  when  performing  system 
identification  is  of  the  utmost  importance.  This  is  especially  true  of  large  space  structures. 
Insulated  from  the  environment  by  the  air-bearing  suspension  system,  the  DTA  was  in  a 
simulated  zero-gravity  environment.  However,  in  earlier  efforts  to  model  the  test  article, 
extraneous  wiring  (not  in  use  at  the  time)  was  attached  to  the  DTA.  These  wires,  adding 
nonlinear  effects,  caused  many  of  the  system’s  modes  to  change  frequency  and  amplitude. 
When  obtaining  an  accurate  model  of  a  system  on  which  one  plans  to  employ  active  control 
strategies,  the  structure  should  be  identified  in  the  same  configuration  in  which  it  will  be 
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utilized.  This  was  the  case  in  this  identification  and  the  removal  of  the  instrumentation  wiring 
proved  valuable  in  obtaining  an  accurate  model  of  the  DTA. 

The  ERA  algorithm  was  validated  first  on  the  truth  model,  a  linear  system.  It  provided 
an  excellent  fit  of  the  truth  model,  even  when  a  15%  noise  was  added  to  the  frequency  response 
functions.  However,  the  same  results  were  not  obtained  from  the  more  complex  PACOSS 
DTA.  There  are  several  reasons  for  this  discrepancy.  First,  the  truth  model  is  a  linear  system, 
whereas  the  DTA  is  a  nonlinear  system.  Since  the  model  being  developed  is  linear,  it  is 
expected  it  would  identify  the  linear  system  better  than  the  nonlinear  system.  Also,  the  noise 
in  the  truth  model  did  not  obscure  the  shape  of  the  FRF  as  sometimes  happens  in  the  DTA. 
Thus,  the  ERA  routine  had  a  much  easier  time  identifying  the  modes  in  the  truth  model. 
Noise  was  added  to  the  truth  model’s  frequency  response  functions  after  they  had  already  been 
obtained.  Noise  has  many  ports  of  entry  into  the  frequency  response  functions  of  the  DTA, 
especially  when  active  Coulomb  damping  compensation  is  being  applied. 

The  ERA  routine  used  for  this  system  identification  was  written  in  MATLAB^^.  A 
Sparc  10  with  dual  55  MHz  processors  was  also  used.  The  Sparc  10  provided  a  needed 
increase  in  speed  and  random  access  memory,  over  the  Sparc  2,  in  identifying  the  8  x  8  MIMO 
system  with  80  states.  The  identification  process  was  limited  to  the  first  40  modes  due  to  the 
limitation  on  computer  memory.  If  a  machine  with  more  memory  had  been  used,  the  number 
of  modes  identified  could  have  been  increased.  If  the  amount  of  memory  available  could  be 
increased,  the  identification  could  be  extend  up  to  20  Hz  at  a  minimum. 

After  obtaining  the  frequency  response  functions,  it  took  more  than  one  run  of  the  ERA 
routine  to  acquire  an  acceptable  model  of  the  DTA.  The  method  first  used  in  determining  the 
number  of  states  to  be  identified  was  to  count  the  number  of  peaks  in  the  FRF.  This  method 
led  to  an  inaccurate  and  erroneous  model.  In  utilizing  this  method,  40  states  were  identified, 
whereas  80  states  should  have  been  identified  below  14  Hz.  This  was  remedied  when  it  was 
discovered  that  there  were  40  modes  (i.e.,  80  first-order  states)  below  the  cutoff  frequency  of 
14  Hz.  This  was  determined  by  looking  at  the  finite  element  model  generated  by  Capt.  Scott 
George  (4).  Learning  from  this  mistake,  it  is  recommended  that  a  modal  analysis  of  some  type 
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be  completed  prior  to  utilizing  the  Eigensystem  Realization  Algorithm.  Doing  so  will  save 
time  and  result  in  a  more  accurate  description  of  the  systems  low  frequency  dynamics. 

As  can  be  seen  in  the  results  of  the  ERA  fit  of  the  DTA,  ERA  fits  the  poles  in  the  system 
rather  well,  but  is  much  more  inaccurate  in  fitting  the  zeros.  Therefore,  it  is  recommended 
that  parameter  optimization  be  applied  to  this  system.  The  A,  B,  C,  and  D  matrices  identified 
by  the  ERA  algorithm  could  be  used  as  starting  points  in  the  time  domain  based  parameter 
optimization  routines  to  improve  the  system  response  fits.  MATLAB’s  System  Identification 
toolbox  contains  some  parameter  optimization  routines  that  could  be  used  to  possibly  obtain 
a  more  accurate  state-space  representation  of  the  PACOSS  DTA’s  low  frequency  dynamics. 
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Appendix  A.  State-Space  Model  FRFs  ofPACOSS  DTA 


This  appendix  is  divided  into  eight  separate  appendices  with  each  appendix  representing 
a  column  of  the  transfer  function  matrix.  These  appendices  contain  the  output  (a  pictorial 
representation)  of  the  frjfilt.m  routine. 

This  appendix  contains  the  eight  measurements  resulting  from  the  exciting  of  Actuator 
1  with  a  random  noise. 
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Figure  A.2  Measurement  2  Excitation  Actuator  1 
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Figure  A.5  Measurement  5  Excitation  Actuator  1 
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Appendix  B.  State-Space  Model  FRFs  ofPACOSS  DTA  (Cont.) 

The  eight  plots  in  this  appendix  are  the  measurements  taken  when  Actuator  2  was 
excited  with  a  random  noise. 
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Figure  B.3  Measurement  3  Excitation  Actuator  2 
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Figure  B.4  Measurement  4  Excitation  Actuator  2 
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Figure  B.6  Measurement  6  Excitation  Actuator  2 
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Appendix  C.  State-Space  Model  FRFs  ofPACOSS  DTA  (Cont.) 

The  eight  plots  in  this  appendix  are  the  measurements  taken  when  Actuator  3  was 
excited  with  a  random  noise. 
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Figure  C.3 
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Figure  C.6  Measurement  6  Excitation  Actuator  3 
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Appendix  D.  State-Space  Model  FRFs  ofPACOSS  DTA  (Cont.) 

The  eight  plots  in  this  appendix  are  the  measurements  taken  when  Actuator  4  was 
excited  with  a  random  noise. 
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Figure  D.3  Measurement  3  Excitation  Actuator  4 
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Figure  D.6  Measurement  6  Excitation  Actuator  4 
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Appendix  E.  State-Space  Model  FRFs  ofPACOSS  DTA  (Cont.) 

The  eight  plots  in  this  appendix  are  the  measurements  taken  when  Actuator  5  was 
excited  with  a  random  noise. 
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Figure  E.3  Measurement  3  Excitation  Actuator  5 
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Figure  E.6  Measurement  6  Excitation  Actuator  5 
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Figure  E.7  Measurement  7  Excitation  Actuator  5 


Appendix  F.  State-Space  Model  FRFs  of  PACOSS  DTA  (Cont.) 

The  eight  plots  in  this  appendix  are  the  measurements  taken  when  Actuator  6  was 
excited  with  a  random  noise. 
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Figure  F.3  Measurement  3  Excitation  Actuator  6 
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Appendix  G.  State-Space  Model  FRFs  ofPACOSS  DTA  ( Cont.) 

The  eight  plots  in  this  appendix  are  the  measurements  taken  when  Actuator  7  was 
excited  with  a  random  noise. 


G-1 


5 


Frequency  (Hz) 


Magnitude 


Figure  G.3  Measurement  3  Excitation  Actuator  7 
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Figure  G.4  Measurement  4  Excitation  Actuator  7 
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Figure  G.6  Measurement  6  Excitation  Actuator  7 
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Figure  G.8  Measurement  8  Excitation  Actuator  7 


Appendix  H.  State-Space  Model  FRFs  of  PACOSS  DTA  (Cont.) 

The  eight  plots  in  this  appendix  are  the  measurements  taken  when  Actuator  8  was 
excited  with  a  random  noise. 
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Figure  H.2  Measurement  2  Excitation  Actuator  8 
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Figure  H.4  Measurement  4  Excitation  Actuator  8 
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Figure  H.5  Measurement  5  Excitation  Actuator  8 
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Appendix  I.  State-Space  Model  FRFs  of  Truth  Model 


The  plots  seen  in  this  appendix  show  the  measured  FRFs  of  the  truth  model,  with 
and  without  noise,  and  their  corresponding  FRFs  generated  by  the  ERA-produced  state-space 
model. 

I.l  Truth  Model  without  Noise 


Figure  1.2  Displacement  Y2  Force  Input  FI 
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Figure  1.4  Displacment  Y 1  Force  Input  F2 
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Figure  1.5  Displacement  Y2  Force  Input  F2 
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Figure  1.8  Displacement  Y2  Force  Input  F3 


Figure  1.9  Displacement  Y3  Force  Input  F3 
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1.2  Truth  Model  with  Noise 


Figure  1. 10  Displacement  Y1  Force  Input  FI 
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Figure  1.13  Displacment  Y 1  Force  Input  F2 
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Figure  1. 14  Displacement  Y2  Force  Input  F2 
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Figure  1.15  Displacement  Y3  Force  Input  F2 


Figure  1. 17  Displacement  Y2  Force  Input  F3 


Figure  1.18  Displacement  Y3  Force  Input  F3 
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Appendix  J.  Computer  Routines 

This  appendix  contains  an  explanation  of  each  of  the  MATLAB^^  ERA  routines  and 
how  they  are  used.  Also,  the  parameter  settings  used  to  identify  the  PACOSS  DTA  and  the 
truth  model  are  presented. 

J.l  ERA  Routines 

The  ERA  routines  used  for  this  identification  can  be  found  in  Appendix  L.  These 
routines,  as  mentioned  earlier,  were  written  by  Capt.  Richard  Cobb.  The  following  sections 
give  a  brief  overview  of  each  routine.  The  last  section  shows  the  settings  used  to  generate  the 
results  shown  in  Appendices  A-H  and  the  corresponding  A,  B,C,  and  D  matrices. 

J.1.1  frfinter.m.  This  function  integrates  the  frequency  response  functions.  The 
data  obtained  by  the  IP  program  was  acceleration  transfer  functions.  When  designing  control 
systems,  acceleration  is  seldom  used;  therefore,  the  frequency  response  functions  must  be 
integrated.  It  follows  then  that  the  number  of  integrations  to  be  performed  depends  on  the 
choice  of  using  either  velocity  or  displacement  FRFs  in  the  control  system  design. 

J.l  .2  mirror. m .  This  function  was  written  in  order  to  produce  a  symmetric  frequency 
response  function  given  an  asymmetric  response  function.  An  FRF  obtained  from  the  DP 
program  (after  ASCII  conversion)  is  input  into  this  m-file.  The  FRF  input  into  ERA  routines  is 
composed  of  (n/2)  +  1  complex  points.  This  points  are  for  the  positive  frequencies.  In  order 
to  perform  an  inverse  fast  Fourier  transform  (IFFT),  an  n  point  symmetric  function  is  needed. 
This  function  mirrors  the  positive  frequencies  into  negative  frequencies;  thus  producing  an 
n  point  symmetric  function.  After  the  mirroring  has  been  accomplished,  the  IFFT  can  be 
performed  on  the  data. 

J.l. 3  mimoera.m.  Given  the  impulse  responses,  number  of  sensors,  number  of 
actuators,  and  number  of  states,  this  function  produces  the  A,  5,  C  and  D  matrices  for  the 
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system.  The  impulse  responses  are  obtained  by  taking  the  inverse  fast  Fourier  transform  of 
the  output  of  mirror.m. 

J.l  A  weave. m  and  mimohank.m.  weave. m  weaves  the  individual  impulse  response 

functions  into  the  proper  form  to  be  utilized  by  mimohank.m. 

mimohank.m  constructs  the  block  Hankel  matrix  for  a  MIMO  system.  This  routine 
needs  the  output  of  weave. m  (the  combined  impulse  functions)  and  the  number  of  sensors  in 
the  system.  The  output  of  this  system  is  the  block  Hankel  matrix  for  use  in  mimoera.m.  A 
singular  value  decomposition  will  be  performed  on  a  portion  of  this  Hankel  matrix. 

J.l. 5  parm.m.  This  routine  is  a  utility  function.  If  pis  the  parameters  vector,  typing 
parin(p)  in  the  MATLAB^^  command  window  produces  output  on  the  screen  showing  the 
current  parameter  settings  and  explains  each  parameter.  By  typing  parm,  the  help  file  for 
parm.m,  which  explains  each  parameter,  is  displayed  on  the  screen. 

J.l  .6  frffilt.m.  This  routine  encompasses  all  of  the  other  routines.  This  function  calls 
mirror.m,  performs  the  inverse  Fourier  transform,  and  (after  some  filtering)  calls  mimoera.m. 
This  routine  needs  the  FRF  matrix,  the  frequency  vector  (FreqV),  and  the  parameters  vector. 
The  output  of  this  function  is  the  A,  5,  C,  and  D  matrices  as  well  as  FRF/ERA  plots. 

J.2  Parameter  Settings 

The  next  two  sections  give  the  parameter  settings  used  in  the  identification  of  the 
PACOSS  DTA  and  the  identification  of  the  truth  model. 

7.2.7  PACOSS  DTA.  The  system  to  be  identified  had  8  inputs  and  8  outputs. 
Looking  to  identify  the  modes  below  14  Hz,  the  number  of  states  was  set  to  80.  The  Ih  factor, 
one  of  the  parameters  which  determines  the  size  of  the  Hankel  matrix,  was  set  to  5  to  lengthen 
the  time  series  used  to  generate  the  Hankel  matrix.  One  integration  was  performed  on  the 
data.  The  original  data  was  acceleration  frequency  response  data.  It  was  determined  that  one 
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integration  to  obtain  velocity  response  data  was  sufficient  to  determine  the  dynamics  of  the 
system.  Two  integrations  could  have  been  performed;  however,  in  the  interest  of  conserving 
CPU  time,  only  one  integration  was  performed.  No  filtering  was  used.  This  was  accomplished 
by  setting  the  low  pass  filter  cutoff  parameter  to  -1.0.  Due  to  the  fact  that  acceleration  data 
was  being  used  and  that  initial  accelerations  are  usually  rather  large,  it  was  found  that  the  first 
three  data  points  needed  to  be  set  to  zero.  Figures  were  plotted  using  a  separate  routine  shown 
in  Appendix  L.  A  discrete  model  was  used.  The  number  of  sample  points  was  set  to  1024. 
It  was  also  decided  that  a  singular  value  cutoff  was  not  needed.  These  parameter  settings 
were  input  into  the  parameters  vector,  which  is  called  by  the  ERA  routines  to  generate  the 
state-space  model  of  the  PACOSS  DTA  found  in  this  document. 

J.2.2  Truth  Model.  The  parameter  settings  used  for  the  truth  model  are  not  much 
different  from  those  used  in  the  identification  of  the  test  article.  Only  those  that  differ  will  be 
discussed  here. 

The  truth  model  system  has  eight  states  with  3  inputs  and  3  outputs.  No  integrations 
were  used  since  the  output  of  the  system  was  displacement.  The  Ih  factor  was  set  to  4  to 
increase  the  size  of  the  Hankel  matrix  since  the  system  was  small.  Also,  no  zeros  were  placed 
at  the  beginning  of  the  data  sets.  The  remaining  settings  were  the  same  as  those  listed  for  the 
PACOSS  DTA. 

In  order  to  generate  the  state-space  model  and  provide  for  user-interface,  a  routine 
(erasetup.m)  was  written  which  sets  up  the  FRF  matrix,  clears  memory  to  speed  processing, 
and  saves  the  workspace.  Another  function  iplotfrfs.m)  was  written  which  loads  the  workspace, 
plots  the  FRFs  and  their  corresponding  ERA  fits,  and  saves  them  automatically  to  a  file.  The 
code  for  both  of  these  routines  is  shown  in  Appendix  L. 
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Appendix  K.  C  Code  Used  in  Friction  Compensation 


K.l  C  Code  for  Coulomb  Damping  Compensator 

This  section  contains  the  code  use  to  implement  the  Coulomb  damping  compensator. 


/***************************************************************** 

* 

*  COULOMB  DAMPING  COMPENSATOR 

* 

*  Author:  Capt  Robert  Daryl  Woods  Date:  29  Aug  1994 

* 

*  This  routine  was  written  to  compensate  for  the  Coulomb 

*  Damping  present  in  each  actuator.  The  Coulomb  damping 

*  coefficients  have  already  been  determined  through  experi- 

*  mentation  and  tweaked  using  the  Optima  to  ensure  the 

*  compensation  voltage  being  applied  does  not  start  the 

*  actuator.  The  program  has  to  be  recompiled  and  downloaded 

*  each  time  the  disturbance  actuator  changes. 

* 

*  USAGE 

* 

*  The  program  is  compiled  using  the  couldamp.cmd  script 

*  file.  This  file  looks  like: 

* 

*  cc  -c  -0  -f 68881  -0  couldamp.o  couldamp.c 

*  Id  -N  -0  couldamp  couldamp.o  convrtl6.o  qc.l  std.l  stdio.l 

* 

*  The  file  is  than  downloaded  to  the  Optima  using: 

•k 

*  download  couldamp 

* 

*****************************************************************/ 


/*  DEFINITIONS  */ 

#  define  M  9 

/*  INCLUDE  STATEMENTS  */ 

#include  "stdio.h" 
tinclude  "math.h” 
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float  *APalloc(); 


/*  MAIN  PROGRAM  */ 

main ( ) 

{ 

float  *b,  *C,  *c_force,  *offset,  f_out; 

float  dt,  zone,  incal, out cal, endtime, endcount, in_gain, out_gain; 
int  i,  ii,  j, k,  1; 

/*  Dimensionalize  Variables  */ 

APInit (0) ; 

b  '  =  APalloc  (M) ; 

C  =  APalloc  (M) ; 

c_force  =  APalloc  (M)  ; 
offset  =  APalloc  (M)  ; 


* 

Set  Offsets  for 

each  Channel 

*/ 

offset [0] 

=  0.051; 

/* 

Offset  for  Channel 

1 

*/ 

offset [1] 

=  0.010; 

/* 

Offset  for  Channel 

2 

*/ 

offset [2] 

=  0.071; 

/* 

Offset  for  Channel 

3 

*/ 

offset [3] 

=  0.071; 

/* 

Offset  for  Channel 

4 

*/ 

offset [4] 

=  0.062; 

/* 

Offset  for  Channel 

5 

*/ 

offset [5] 

=  0.064; 

/* 

Offset  for  Channel 

6 

*/ 

offset [6] 

=  0.049; 

/* 

Offset  for  Channel 

7 

*/ 

offset [7] 

=  0.046; 

/* 

Offset  for  Channel 

8 

*/ 

offset [8] 

=  0.036; 

/* 

Offset  for  Channel 

9 

*/ 

•k 

Set  Coulomb  Damping  Amounts 

*/ 

c  force [0] 

=  0.0325 

/* 

Actuator  1  LVT 

*/ 

c_force [1] 

=  0.0500 

/* 

Actuator  2  LVT 

*/ 

c_force [2] 

=  0.0250 

/* 

Actuator  3  LVT 

*/ 

c  force [3] 

=  0.0100 

/* 

Actuator  4  LVT 

*/ 

c  force [4] 

=  0.0150 

/* 

Actuator  5  LVT 

*/ 

c  force [5] 

=  0.0500 

/* 

Actuator  6  LVT 

*/ 

c  force [6] 

=  0.0325 

/* 

Actuator  9  LVT 

*/ 

c  force [7] 

=  0.0175 

/* 

Actuator  10  LVT 

*/ 

c  force [8] 

=  0.00; 

* 

Assign  Conversions  and  Variables 

*/ 

incal  =  1./0.6;  /*  Convert  input  Volts  ->  in/sec  */ 
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outcal  =  1./0.2;  /*  Convert  output  lbs  ->  Volts  */ 

dt  =  0.005;  /*  Sample  Time  */ 

zone  =  0.2;  /*  Dead  Zone  level  (+/-  in/sec)  */ 

endtime  =  3600;  /*  End  of  Data  Taking  */ 

endcount  =  endtime/dt;  /*  End  of  Loop  */ 

in_gain  =  -0.0003;  /*  Gain  on  input  (ATD)  Board  */ 

out_gain  =  3333.333;  /*  Gain  on  Output  (DTA)  Board  */ 

/* - * 

*  -  Intialize  cards  with  dt  (second)  sample  time 

*  - *  j 

if  (CVT_INIT(dt,3) )  { 
printf ("Sample  Time  Out  of  Range\n"); 
exit (0) ; 

} 

/* - * 

*  -  Start  Conversion  Clock 

*  - */ 

CVT_STRT ( ) ; 

I  * - * 

*  -  Start  For  Loop 

*  - - *  j 

for  (i=0;i  <=  endcount;  i++)  {  /*  Start  loop  */ 

I  * - * 

*  -  Get  Analog  Input 

*  - */ 

if  (CVT_ATD(&C[0] ,  &C[1],  &C[2],  &C[3],  /*  FirstTCard  Data  */ 

&C[4],  &C[5],  &C[6],  &C[7],  /*Second  Card  Data  */ 
&C[8],  &C[9] ,&C[10],  &C[11]))  {  /*  Third  Card  Data  */ 
printf ("Overrun  Error  \n") ; 
exit (0)  ; 

} 

/* - * 

*  -  Multiply  input  by  Input  gain  &  subtract  offset 

*  - */ 

for  (j=0;  j<=8;  j++)  { 

C[j]  =  (C[j]  *  in_gain)  -  offset  [j]; 

} 

/* - * 

*  -  Convert  1st  8  inputs  to  Velocity 
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* 


*/ 

for  (1=0;  1<=7;  1++)  { 

C[l]  =  C[l]  *  incal; 

} 

/* - * 

*  -  SIGNUM  CALCULATIONS 

*  - */ 

for  (ii=0;  ii<=7;  ii++)  { 

f_out  =  0 . ; 

if  (C[ii]  >  zone)  f_out  =  c_force[ii]; 
if  (C[ii]  <  -zone)  f  out  =  -c  force [ii]; 


b[ii]  =  f_out  *  0.9;  /*  90%  of  friction  force  */ 

} 

/* - * 

*  -  Convert  lbs  to  Voltage 

*  - - - : - */ 

for  ()c=0;  k<=7;  ]c++)  { 
b[k]  =  b[k]  *  outcal  *  out_gain; 

} 

j  -k - k 


*  -  Add  compensation  to  random  input  from  Analyzer 

*  - */ 

b[8]  =  (C[8]  *  out_gain)  +b[0]; 

j  k - * 

*  -  Send  Analog  Output 

*  - *  j 

if  (CVT_DTA(&b[0] ,  &b[l],  &b[2],  &b[3],  /*  First  Card  Data  */ 

&b[4],  &b[5],  &b[6],  &b[7],  /*Second  Card  Data  */ 
&b[8],  &b[9],&b[10],  &b[ll]))  {  /*  Third  Card  Data  */ 
print f ("Overrun  Error  \n")  ; 
exit  (0)  ; 

} 

}  /*  End  Loop  */ 

}  /*  End  Program  */ 
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K.2  Coulomb  Damping  “Tweaking”  Routine 


This  routine  was  used  to  “tweak”  the  calculated  Coulomb  damping  coefficients.  The 
compensated  output  is  displayed  on  the  IP  program  to  determine  the  accuracy  of  the  friction 
coefficient  for  each  actuator.  Depending  on  the  output,  the  coefficient  for  that  actuator  is 
adjusted  either  up  or  down  and  the  routine  is  re-compiled,  downloaded,  and  executed  again 
on  the  Optima/3.  This  cycle  continues  until  all  of  the  coefficients  are  correct.  This  may  have 
to  be  accomplished  just  prior  to  taking  data  for  a  compensated  system  since  the  coefficients 
seem  to  be  affected  by  the  temperature  and  humidity  changes. 

* 

*  Friction  Tweaking  Routine 

"k 

*  Author:  Capt  Robert  Daryl  Woods  Date:  23  Aug  1994 

* 

*  This  routine  was  written  to  tweak  the  dry  friction 

*  constants  on  Actuators  9  &  10.  It  can  also  be  used  on  the 

*  other  6  actuators.  The  LVT  output  will  be  looked  at  to 

*  determine  the  friction  constant.  The  free  decay  method 

*  (lifting  mass  to  a  preset  height  and  releasing)  will  be 

*  employed.  The  compensated  output  will  be  looked  at  on  the 

*  IP  display  to  determined  the  accuracy  of  the  coefficient. 

*  The  coefficient  is  then  adjusted  up  or  down  as  appropriate. 

* 

kkkk-kkk-k-k-k'k'kk'k-k'kkkk-kk'kkk’kk'kk'kkk'k'kic'kk’kk'k'k'k'kk'kkkk'kk'kkk'kickkkkic'kk'k'k'k'k^ 

/*  DEFINITIONS  */ 

#  define  M  9 

/*  INCLUDE  STATEMENTS  */ 

linclude  "stdio.h" 

#include  "math.h" 

float  *APalloc(); 

/*  MAIN  PROGRAM  */ 

main ( ) 

{ 

float  *b,  *C,  *c_force,  *offset,  f_out; 
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float  dt, zone, incal, outcal, endtime, endcount, in_gain, out_gain; 


int  i,ii, j,k,l; 

/*  Dimensionalize  Variables  */ 

APInit(O) ; 

b  =  APalloc  (M)  ; 

C  =  APalloc  (M)  ; 

c_force  =  APalloc  (M) ; 
offset  =  APalloc  (M)  ; 

*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 


offset [0] 
offset [1] 
offset [2] 
offset [3] 
offset [4] 
offset [5] 
offset [6] 
offset [7] 
offset [8] 


0.051 

0.010 

0.071 

0.071 

0.062 

0.064 

0.049 

0.046 

0.036 


Set  Offsets  for  each  Channel 

/*  Offset  for  Channel  1 
/*  Offset  for  Channel  2 
/*  Offset  for  Channel  3 
/*  Offset  for  Channel  4 
/*  Offset  for  Channel  5 
/*  Offset  for  Channel  6 
/*  Offset  for  Channel  7 
/*  Offset  for  Channel  8 
/*  Offset  for  Channel  9 


Initial  Coulomb  Damping  Coefficients 


c_force [0] 
c_force [1] 
c_force [2] 
c_force [3] 
c_force [4] 
c_force [5] 
c_force [6] 
c  force [7] 


0.0591 

0.0565 

0.0549 

0.0353 

0.0574 

0.0518 

0.2824 

0.3188 


Actuator  1  LVT 
Actuator  2  LVT 
Actuator  3  LVT 
Actuator  4  LVT 
Actuator  5  LVT 
Actuator  6  LVT 
Actuator  9  LVT 
Actuator  10  LVT 


*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 


/*  Assign  Conversions  and  Variables  */ 

incal  =  1./0.6;  /*  Convert  input  Volts  ->  in/sec  */ 

outcal  =  1./0.2;  /*  Convert  output  lbs  ->  Volts  */ 

dt  =  0.005;  /*  Sample  Time  */ 

zone  =  0.2;  /*  Dead  Zone  level  (+/-  in/sec)  */ 

endtime  =  3600;  /*  End  of  Data  Taking  */ 

endcount  =  endtime/dt;  /*  End  of  Loop  */ 

in_gain  =  -0.0003;  /*  Gain  on  input  (ATD)  Board  */ 

out_gain  =  3333.333;  /*  Gain  on  Output  (DTA)  Board  */ 
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/* - : - * 

*  -  Intialize  cards  with  dt  (second)  sample  time 

*  - */ 

if  (CVT_INIT{dt,3) )  { 

printf ("Sample  Time  Out  of  Range\n"); 
exit  (0) ; 

} 

/* - * 

*  -  Start  Conversion  Clock 

*  - - - */ 

CVT  STRT ( ) ; 


/* - 

*  _ 

* _ 

for  (i=0;i 


Start  For  Loop 


<=  endcount;  i++)  { 


_ * 

- */ 

/*  Start  loop  */ 


/* 


*  -  Get  Analog  Input 

*  - */ 

if  (CVT_ATD(&C[0],  &C[1],  &C[2],  &C[3],  /*  First  Card  Data  */ 

&C[4],  &C[5],  &C[6],  &C[7],  {  /*Second  Card  Data  */ 

printf ("Overrun  Error  \n")/ 
exit (0) ; 

} 

/* - * 

*  -  Multiply  input  by  Input  gain  &  subtract  offset 

*  - */ 


for  (j=0;  j<=7;  j++)  { 

C[j]  =  (C[j]  *  in_gain)  -  offset [j]; 

} 

/* - * 

*  -  Convert  1st  8  inputs  to  Velocity 

*  - */ 


for  (1=0;  1<=7;  1++)  { 
C[l]  =  C[l]  *  incal; 


} 

/* - * 

*  -  SIGNUM  CALCULATIONS 

*  - */ 


for  (ii=0;  ii<=7;  ii++)  { 
f  out  =  0 . ; 
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if  (C[ii]  >  zone)  f_out  =  c_force [ii]  ; 
if  (C[ii]  <  -zone)  f_out  =  -c_force [ii] ; 
b[ii]  =  f_out; 

} 

/* - * 

*  -  Convert  lbs  to  Voltage 

*  - */ 

for  (k=0;  k<=7;  k++)  { 

b[k]  =  b[k]  *  outcal  *  out_gain; 

} 

/* - * 

*  -  Send  Analog  Output 

*  - */ 

if  (CVT_DTA(&b[0] ,  &b[l],  &b[2],  &b[3],  /*  First  Card  Data  */ 

&b[4],  &b[5],  &b[6],  &b[7]))  {  /*Second  Card  Data  */ 

print f ("Overrun  Error  \n")  / 
exit (0) ; 

} 

}  /*  End  Loop  */ 

}  /*  End  Program  */ 
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Appendix  L.  MATLAB^^  Code 


This  appendix  is  composed  of  three  sections.  The  first  section  contains  the  m-files 
which  compose  the  ERA  algorithm.  The  next  section  contains  m-files  which  load  the  data  and 
implement  Capt.  Cobb’s  routines  with  the  proper  parameters.  The  last  section  contains  the 
interactive  routine  used  to  calculate  the  Coulomb  damping  coefficients. 

L.1  MATLAB™  ERA  Code 

This  section  contains  m-files  written  to  perform  the  Eigensystem  Realization  Algorithm. 
The  first  section  contains  the  actual  ERA  routine.  These  files  were  written  by  Capt.  Richard 
Cobb,  a  Ph.D.  candidate  at  the  Air  Force  Institute  of  Technology.  A  description  of  each  m-file 
is  given  in  Appendix  J. 

L.l  .1  frffiltm. 


function  [a,b,c,d,Ts]=frffilt  (FRF,FreqV,paramtrs) 

% 

%  [a,b,c,d, Ts]=filtfrf (FRF,FreqV,pararatrs) 

%  FRF=[frfll  frf21  ...  frfl2  frf22 . ] 

%  frfij  =  frequency  response  function  at  sensor  i  from  input  j 
%  FreqV  =  frequency  vector  in  Hz,  starting  at  0  Hz . 

%  paramtrs  =  parameter  vector  as  defined  below 
% 

%  paramtrs= [nmeas, ninput , lh,ns, ni, ford, Hzbreak, nz, fn, cord, cut] 

% 

%  nmeas  =  #  of  sensors 
%  ninput  =  #  of  actuators 

%  Ih  is  a  factor  2  <  Ih  <  length  (FRF) /ns  to  determine  size  of 
%  Hankel  matrix  typically  3-8.  The  resulting  Hankel  matrix 

%  is  nxn  where  n  is  approximately  ns*nmeas*ninput*lh/2 

%  ns  =  #  number  of  states  for  model 

%  ni  =  #  of  integrations  on  data  (set  to  1  for  accel  to  velocity) 

%  ford  =  butterworth  filter  order,  ignored  if  break  is  negative 

%  Hzbreak  =  cutoff  for  low  pass  filter  of  impulse  response  in  Hz 
%  set  break  <  0  for  no  filtering 

%  nz  =  #  of  zeros  at  low  freq  end  of  spectrum 
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% 


fn  =  figure  number  for  first  plot  (<0  to  surpress  plots) 

%  cord  =  1  for  continuous  model  or  2  for  discrete  model 
%  cut  =  1  to  prompt  for  singular  value  cutoff  (using  figure  fn) 

%  ndpts  =  number  of  original  sample  points 

%  (ie  2"N,  512,1024,2048)  set  to  zero  for  no  padding 

% 

%  rcobb  15/feb/94 

lenfrf=size (FRF, 1) ; 
if  length (FreqV)  -  lenfrf  ~=  0 

error (' incompatable  length  FRF  and  Freq  vectors') 
end 

nmeas=paramtrs (1) ;  ninput=paramtrs  (2) ;  lh=paramtrs (3) ; 
ns=paramtrs (4) ;  ni=paramtrs (5) ;  ford=paramtrs (6) ; 
Hzbreak=paramtrs (7) ;  nz=paramtrs (8) ;  fn=paramtrs ( 9) ; 
cord=paramtrs (10) ;  cut=paramtrs (11) ;  ndpts=paramtrs (12) ; 

%  Integrate  frf  to  get  vel/force  or  disp/force  xfer 
%  function 
for  i=l :ni 

FRF=frfinter(FRF,FreqV*2*pi)  ; 
end 

%  pad  frf  with  zeros  to  match  sample  rate 
if  ndpts  '=  0 

FRF= [FRF; zeros (ndpts+l-lenfrf , size (FRF, 2) ) ] ; 
end 

%  compute  sample  rate  of  data 

Ts=l/( (FreqV(2)-FreqV(l)) *(size(FRF,l)-l)*2) ; 

%  set  low  freq  end  to  zero 
if  nz  >  0 

for  i=l ;nmeas*ninput 
FRF (1 :nz, i) =zeros (nz, 1) ; 
end 
end 

%comput  impulse  response  and  filter  data 
for  i=l :nmeas*ninput 
G( : , i) =mirror (FRF ( : , i) ) ; 
end 

h=real(ifft(G) ) ; 
if  Hzbreak  >  0 

[B, A] =butter ( ford, Hzbreak*2*Ts ) ; 
for  i=l :nmeas*ninput 

h  ( : ,  i)  =filtfilt  (B,  A,h( : ,  i) )  ; 
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end 

end 

% 

[a,b,c,d]=mimoera(h(2:ns*lh+l, :) ,nmeas,ninput,ns,fn,cut) ; 

Q. 

“O 

if  cord  ==  1 

[a,b]=d2c(a,b,Ts) ;  Ts=[]; 

%  cut  identified  poles  and  zeros  above  max  frequency 

%  [a,b, c, d] =cutzeros (a,b, c, d,max (FreqV) )  ; 

end 

if  fn  >  0  %  print  plots  for  positive  figure  number 

for  i=l:ninput 
if  cord  ==  1 

[mag,phase] =bode (a,b, c, d, i,FreqV*2*pi) ; 
else 

[mag, phase ] =undbode (a,b, c, d,Ts, i,FreqV*2*pi) ; 
end 

figure  (fn-l+i) ,  subplot  (2,1,1),  ... 
semilogy(FreqV,abs(FRF(l:lenfrf,nmeas*(i-l)+l:  .  . . 

nmeas* (i-1) +nmeas) ) , FreqV, mag) 
figure  (fn-l+i) ,  subplot  (2,1,2),  ... 
plot (FreqV, (angle (FRF (1 :lenfrf,nmeas* (i-1) +1 :  ... 
nmeas* (i-1) +nmeas) ) ) *180/pi, FreqV, phase) 

end 

end 

%eof 
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L.1.2  mimohank.m. 


function  H=mimohank (Y, nmeasure) 

% 

%  H=mimohank (Y,  #  of  sensors) 

% 

%  This  routine  constructs  the  block  Hankel  matrix  for  a 
%  MIMO  system.  Y  is  the  output  of  WEAVE  (combined  impulse 
%  functions) 

%  rcobb  3/feb/94 

% 

nin=nmeasure; 

%  correct  length  of  Y  to  avoid  diagonal  conflicts 
[x, ncol] =size (Y)  ; 

Y=Y  (1  :fix  (x/  (4*nin)  )  *nin*4,  : )  / 

[x,ncol]=size(Y)  ; 
nrows=x/nin; 
nrblock=nrows/2; 
ncblock=nrows/2; 

H=zeros (nrblock*nin, ncblock*ncol) ; 
for  j=0 :ncblock-l 

H(l:nrblock*nin,ncol* j+l:ncol*j+ncol)=Y( j*nin+l:  . .  . 

j*nin+nin*nrblock, 1 :ncol) ; 

end 

%eof 
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L.13  mimoera.m. 


function  [ amod, bmod,  cmod,  dmod]  =minioera  (h,  nmeas ,  ninput ,  nstates,  f n,  cut ) 
% 

%  [a,b,c,d]=mimoera(  [hi  h2  of  sensors,#  of  actuators, 

%  #  of  states, cut) 

%  where 

%  hi  =  impulse  response  to  input  i 

%  hi  has  as  many  columns  as  there  are  sensors 

% 

%  fn,cut  =  if  cut  =  1  then  a  singular  value  plot  will  be 

%  displayed  in  figure  (fn)  and  prompt  the  user 

% 

%  rcobb  15/feb/94 


if  nmeas *ninput  ==  1 

Y=h; 

else 

Y=weave (h, nmeas, ninput) ; 
end 

H=mimohank (Y, nmeas) ; 
n=size (H, 1) /nmeas; 

Ho=H (1 : (n-1) *nmeas, 1 : (n-1) *ninput) ; 

H1=H (1 : (n-1) * nmeas, ninput +1 :n*ninput) ; 

disp (sprintf (' SVD  of  a  %6i  by  %6i  matrix  in  progress' , size (Ho) ) ) 
[U,S,V]=svd(Ho,0) ; 

% [U2,psi,V2]=svdra (Ho, nstates,  3) ; 

% [U4, S, V4]=svd(psi) ; 

%U=U2*U4; 

%V=V2*V4; 

s=diag(S) ; 

if  length (s)  <  nstates 

error ('increase  length  of  time  series') 

elseif  length (s)  <  2*nstates 

disp (' Warning,  increase  length  of  time  series') 
end 

%  plot  singular  values  and  prompt  for  cutoff 
if  cut  ==1 
figure  (fn) 
clf;  semilogy(s) 
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xlabel (' Select  the  Singular  Value  Cutoff  (#  of  states)') 
ylabel('log  of  sv' ) 

nstates=round(myginput) ;nstates=nstates (1) ; 
disp (sprintf (' %5i  states  selected' ,nstates) ) 
end 

s=s (1 :nstates) ; 
u=U ( : , 1  instates) ; 
v=V ( : , 1  instates) ; 

Ep= [eye (ninput)  zeros (ninput, (n-2) *ninput) ] ' ; 

Eq=  [eye  (nmeas)  zeros  (nmeas,  (n-2)  *nineas)  ] ; 
srt=sqrt (s) ; 

amod=diag ( 1 . / srt) *u' *Hl*v*diag (1 . / srt) ; 
bmod=diag(srt) *v' *Ep; 
cmod=Eq*u*diag(srt)  ; 
dinod=zeros  (nmeas, ninput) ; 

%eof 
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L.1.4  frfinter.in. 


function  [frfv]=frfinter  (frf,w) 

% 

%  [frfv]=frfinter  (frf,w) 

% 

%  This  function,  which  integrates  an  frf=[frfl  frf2  ...], 

%  is  used  to  convert  accelerometer  (or  velocity)  measurements 
%  to  velocity  (or  position)  measurements.  W  is  the  frequency 
%  vector  in  rad/sec. 

%  rcobb  2/feb/94 

[m,n]=size(frf) ; 
for  k=l :n 

for  i=2 :m 

frfv (i, k) =frf (i, k) / (sqrt (-1) *w(i) ) ; 

end 

end 

for  k=l:n 
if  w(l) ==0, 
frfv(l,k) =0; 
else 

frfv(l,k) =frf (l,k) / (sqrt (-1) *w(i) ) ; 

end 

end 

%eof 
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L.1.5  mirror.m. 


function  X=mirror(Y) 


X=mirror (Y) 


%  This  function  produces  an  'n'  point  symmetric  FFT  from 

%  an  (n/2+1)  point  asymmetric  FFT. 

%  R.  Cobb,  l/Feb/94 


n=size (Y, 1) ; 


for  i=l:n-2 

Ym(n-i)  =  Y (i+1) ; 

end 


X=[Y;Ym(2:n-l) '  ]  ; 
%  eof 
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L.1.6  weave.in. 


function  Y=weave(X,nmeas,nin) 

% 

%  Y=weave ( [hll  h21...  hl2  h22..],#  of  sensors,#  of  actuators) 
% 

%  This  function  weaves  the  impulse  response  vectors  together 
%  for  use  by  the  MIMO  Hankel  algorithm  where 
% 

%  h(ij)  =  impulse  response  to  input  j  from  sensor  i 

% 

%  see  MIMOHANK  ,MIMOERA 

%  rcobb  3/feb/94 

nr=size (X, 1) ; 

Y=zeros (nr*nmeas, nin) ; 
for  i=0 :nr-l 
for  k=0 :nin-l 
for  j=0:nmeas-l 

Y(i*nmeas+l+j,k+l) =X(i+l,k*nmeas+j+l) ; 
end 
end 
end 
%eof 
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L.1.7  parin.in. 


function  parm(p) 

%  paramtrs= [nmeas, ninput, Ih, ns,ni, ford, Hzbreak, nz,  fn,  ... 

%  cord, cut, ndpts] 

% 

%  nmeas  =  #  of  sensors 
%  ninput  =  #  of  actuators 

%  Ih  is  a  factor  2  <  Ih  <  length  (FRF) /ns  to  determine  size  of 
%  of  Hankel  matrix  typically  3-8.  The  resulting  Hankel 

%  matrix  is  nxn  where  n  “=  ns*nmeas*ninput*lh/2 

%  ns  =  #  number  of  states  for  model 

%  ni  =  #  of  integrations  on  data  (set  to  1  for  accel  to  vel.) 

%  ford  =  butterworth  filter  order,  ignored  if  break  is  negative 

%  Hzbreak  =  cutoff  for  low  pass  filter  of  impulse  response 
%  in  Hz  (set  break  <  0  for  no  filtering) 

%  nz  =  #  of  zeros  at  low  freq  end  of  spectrum 
%  fn  =  figure  number  for  first  plot  (<0  to  surpress  plots) 

%  cord  =  1  for  continuous  model  or  2  for  discrete  model 

%  cut  =  1  to  prompt  for  singular  value  cutoff  (using  Fig.  fn) 

%  ndpts  =  number  of  original  sample  points 

%  (ie  2''N,  512,1024,2048)  set  to  zero  for  no  padding 

if  nargin  ==  0 
help  parm 
else 

disp (sprintf (' %5i  =  #  of  sensors' ,p (1) ) ) 
disp (sprintf (' %5i  =  #  of  actuators' ,p (2) ) ) 

disp  (sprintf  (' %5i  =  Ih  is  a  factor  2  <  Ih  <  length  (FRF) /ns  to  .  . . 

determine  size  of  hankel  matrix  ',p(3))) 
disp (sprintf (' %5i  =  #  number  of  states  for  model' ,p (4) ) ) 
disp  (sprintf  (' %5i  =  #  of  integrations  on  data  (set  to  1  for  ... 
accel  to  velocity) ' ,p (5) ) ) 

disp  (sprintf  (' %5i  =  butterworth  filter  order,  ignored  if  break  .  .  . 

is  negative' ,p (6) ) ) 

disp  (sprintf  (' %5 .  If  =  cutoff  for  low  pass  filter  of  impulse  ... 

response  in  Hz  set  break  <  0  for  no  .  . . 
filtering' ,p  (7) ) ) 

disp (sprintf (' %5i  =  #  of  zeros  at  low  freq  end  of  spectrum' ,p (8) ) ) 
disp  (sprintf  (' %5i  =  figure  number  for  first  plot  (<0  to  . .  . 

suppress  plots) ' ,p (9) ) ) 

disp  (sprintf  (' %5i  =  1  for  continuous  model  or  2  for  discrete  ... 
model' ,p (10) ) ) 
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disp (sprintf ( ' %5i  =  1  to  prompt  for  singular  value  cutoff  .  .  . 

(using  figure  fn)',p(ll))) 

disp (sprintf (' %5i  =  number  of  original  sample  points  set  to  .  .  . 

zero  for  no  padding' ,p (12) ) ) 

end 

%eof 
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L.2  MATLAB^^  ERA  Setup  Routine 


The  routines  presented  in  this  section  were  written  to  setup  the  FRF  vector,  to  clear 
variables  and  vectors  not  needed  for  computation,  and  to  start  fiffilt.m  to  obtain  the  state  space 
model. 

L.2.1  erasetup.m. 


% - 

%  ERASETUP 


%  This  routine  ERASETUP  runs  the  ERA  algorithm.  It  allows 

%  (prompts)  the  user  to  input  the  individual  parameters  which 
%  make  up  the  parameters  vector  (p)  needed  by  FRFFILT.  It 
%  assembles  the  parameters  into  vector  p,  loads  the  transfer 
%  functions  (xfer**'s)  using  FRFBUILD  and  starts  the  ERA  routine 
%  in  the  FRFFILT  function. 

%  Daryl  Woods 

% - - 

% 

%  Prompt  user  to  pick  parameters  % 

% 

nmeas  =  input  ('Number  of  measurements:  '); 
ninput  =  input ('Number  of  actuators:  '); 

Ih  =  input ('Input  Ih  (determines  size  of  Hankel  matrix) :  '); 
ns  =  input  ('Number  of  states  for  the  model:  '); 
ni  =  input ( 'Number  of  intergrations  on  the  data:  '); 
ford  =  input  ('Order  of  Butterworth  Filter:  ' )  ; 

Hzbrk  =  input (' Cutoff  for  the  Low-Pass  Filter  ... 

(<0  for  no  filtering)  :  ' ) ; 

nz  =  input  (' Number  of  zeros  at  low  frequency  end  of  .  .  . 
spectrum:  ' )  ; 

fn  =  input  ('Figure  number  for  first  plot  (<0  to  suppress  .  .  . 
plots) :  ' ) ; 

cord  =  input (' 1-Continuous  or  2-Discrete:  '); 
cut  =  input  ('Set  equal  to  1  for  SV  cutoff:  ' ) ; 
ndpts  =  input ('Number  of  original  points:  ' )  ; 

% 

%  Setup  parameters  vector  % 

p  =  [nmeas, ninput, Ih, ns, ni, ford, Hzbrk, nz, fn, cord, cut, ndpts] ; 

% 

%  Prompt  user  to  see  the  Parameters  vector  % 
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disp(' Would  you  like  to  see  the  parameters  vector  (p) ' ) ; 
yorn  =  input ('1  for  Yes/2  for  No  '); 

% 

if  yorn  ==  1; 
parm(p) 

end 

% 

%  Clear  individual  parameters  since  they  are  in  p  vector  % 
clear  yorn  nmeas  ninput  Ih  ns  ni  ford  Hzbrk  nz  fn  cord; 
clear  cut  ndpts; 

% 

%  Load  in  workspace  with  FRFs  and  build  FRF  matrix  % 
load  xfers3;  %  Loads  FRFs  and  FreqV  % 

frfbuild;  %  Place  FRFs  into  FRF  matrix  % 

clearxfer;  %  Clears  individual  FRFs  to  save  memory  % 

% 

%  Call  frffilt  to  start  ERA  routine% 

[a,b,c,d,  ts]  =  frffilt  (FRF, FreqV, p) ; 

% 

%  End-Of-File 
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L.2.2  frfbuild.m. 


% - 

%  FRFBUILD 


%  This  routine  builds  the  frequency  response  function 
%  (FRF)  matrix  for  use  in  the  FRFFILT  function. 

%  Daryl  Woods 


FRF=[]; 

for  j=l:p(2)  %  Input  index  % 

for  i=l:p(l)  %  Measurement  index  % 

FRF  =  eval ( [' [FRF, xfer' , int2str (i) , int2str ( j) ; 

end 

end 
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L.2.3  clearxfers. 


% 


% 

% 

% 

% 

% 

% 

% 


CLEARXFERS 


Since  the  frequency  response  functions  have  been 
placed  in  the  FRF  matrix,  the  FRFs  will  be  deleted  to 
save  on  memory.  At  the  end  of  the  program,  1  and  k 
(the  counters)  are  also  cleared. 

Daryl  Woods 


for  k  =  p (2)  ; 
for  1  =  p (1)  ; 

eval ( [ ' clear  Xfer' , int2str (1) , int2str (k) ] ) ; 
end 
end 
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L.2.4  plotfrfs.m. 


% - 

%  PLOTFRFS 


% 

%  This  routine  plots  each  FRF  and  its 

%  corresponding  ERA  fit  on  a  separate  plot  and  saves 
%  these  plots  to  an  EPS  file  for  use  with  LaTeX, 

%  Daryl  Woods 


lenfrf  =  length  (FRF)  ;  %  Find  length  of  FRF  % 

FRFv  =  frfinter  (FRF,FreqV*2*pi) ;  %  Integrate  FRF  matrix  % 

k  =  8;  %  k  =  8  for  DTA  % 

%k  =  3;  %  k  =  3  for  truth  model  % 


for  i  =  1 :k; 

[mag, phase]  =  dbode(a,b,c,d, ts,i,FreqV*2*pi) ; 
for  j  =  l:k; 

figure  ( j+(k*  (i-1) )) ,  ... 

subplot (211)  ,semilogy(FreqV(l : lenfrf) ,  abs(FRFv(l:  ... 
lenfrf) ,  j+  (k*  (i-1) ) ) ) , '  -SFreqVd :  . . . 
lenfrf) ,mag(l: lenfrf, j) ' ) ,  ... 
ylabel ('Magnitude' ) ,  ... 
xlabel ('Frequency  (Hz)'),  ... 
legend('-' , 'FRF' 'ERA  Fit' ) ,  ... 
subplot (212) ,plot (FreqV(l : lenfrf) , unwrap (angle (FRFv  .  .  . 

(1 : lenfrf, j+(k* (i-1)))) )*180/pi,'-',  ... 
FreqV(l: lenfrf) ,phase(l:lenfrf, j) ,  ... 

' — '), ylabel ('Degrees' ) ,  ... 
xlabel ('Frequency  (Hz)'); 
orient  tall; 

eval(  ['print  -deps  /rwoods/latex/figs/figure' ,  ... 
int2str ( j) , int2str (i) , ' .eps' ] ) 


end 

end 
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L3  Friction  Coefficient  Calculation  Routine 

% - - - 

%  FC.M 

% - 

%  This  M-File  plots  the  free-decay  of  the  actuator  and 
%  calculates  the  friction  coefficient  for  the  actuator, 

%  Data  was  obtained  via  the  Instrument  Program  and  con- 
%  verted  to  Matlab  format  using  PCMATLNK,  The  two  vectors 
%  used  are  2048x1;  however, not  all  of  this  data  is  needed. 
%  Therefore,  a  lower  and  upper  end  of  usable  data  is 
%  specified.  The  user  is  prompted  to  pick  to  peaks  over 
%  which  the  amplitude  change  is  calculated.  The  Coulomb 
%  friction  coefficient  is  then  calculated  and  output  to 
%  the  screen. 

%  Daryl  Woods 

% - 

min=200;  %  Lower  End  of  Range 

max=400;  %  Upper  End  of  Range 

% 

%  Plot  Time  vs  Data  for  that  range 

plot (TimeV (min:max) , Timei3 (min:max) ./ (0. 6) ) 

% 

%  Prompt  user  for  number  of  points  to  be  fit 
n  =  input  ('Number  of  Points  to  be  Fit?'); 
cycles  =  (n  -  1) ; 

[x,y]  =  ginput(n);  %  User  must  click  on  points 
% 

deltaY  =  y(l)-y(n);  %  Get  delta  Y 

% 

%  Calculate  rps  (frequency  in  radians/sec) 
deltaX  =  x(n)  -  x(l); 
freq  =  cycles  /  deltaX; 
rps  =  2  *  pi  *  freq; 

% 

k  =  1.0000;  %  Spring  Constant  (Ib/in) 

% 

FC  =  (k  *  deltaY)  /  (4  *  cycles  *  rps)  ; 
disp ( ' Coulomb  Friction  is : ' ) , FC 
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